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Abstract 


I  present  a  numerical  technique  to  solve  the  time  independent  Boltzmann 
Transport  Equation  for  the  transport  of  neutrons  and  photons.  The  technique 
efficiently  solves  the  discrete  ordinates  equations  with  a  new  iteration  scheme.  I 
call  this  new  scheme  the  angle  space  distribution  iteration  method  because  it 
combines  a  non-linear,  high  angular-resolution  flux  approximation  within 
individual  spatial  cells  with  a  coarse  angular-resolution  flux  approximation  that 
couples  all  cells  in  a  spatial  mesh.  This  is  shown  to  be  an  efficient  alternative  to 
source  iteration. 

The  new  method  is  implemented  using  the  step  characteristic  and 
exponential  characteristic  spatial  quadrature  schemes.  The  latter  was  introduced 
in  1993  and  has  been  shown  to  be  accurate  for  both  optically  thin  and  optically 
thick  spatial  meshes  and  to  produce  strictly  positive  angular  fluxes. 

The  discrete  ordinates  equations  can  be  solved  using  the  conventional 
source  iteration  method.  However,  it  is  well  known  that  this  method  converges 
prohibitively  slowly  for  optically-thick  problems  with  regions  that  are  dominated 
by  scattering  rather  than  absorption.  The  new  scheme  converges  rapidly  even  for 
such  problems.  Numerical  results  show  that  the  new  scheme  is  reliably  accurate 
for  the  problems  intended,  and  that  it  is  fast  and  efficient  in  use  of  memory. 

The  angle  space  distribution  iteration  method  is  demonstrated  in  slab 
geometry,  for  a  single  energy  group,  using  isotropic  cross  sections,  with 
exponential  and  step  characteristic  spatial  quadratures. 
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A  RAPIDLY-CONVERGING  ALTERNATIVE  TO  SOURCE 
ITERATION  FOR  SOLVING  THE  DISCRETE  ORDINATES 
RADIATION  TRANSPORT  EQUATIONS  IN  SLAB  GEOMETRY 


1.  Introduction 

The  conventional  practice  for  evalnating  the  time  independent  discretized, 
Boltzmann  transport  eqnation  is  the  discrete  ordinates  angnlar  qnadratnre 
method  with  trnncated  Legendre  expansions  representing  the  cross  sections.  In 
discrete-ordinates  S„  approximations  of  large  transport  problems  the  nnderlying 
linear  Boltzmann  problem  is  discretized  in  space  and  angle  and  the  resnlting 
system  of  algebraic  equations  is  solved  iteratively  using  source  iteration.  If  the 
physical  system  contains  regions  that  are  diffusive  and  optically  thick,  source 
iteration  can  be  so  slow  to  converge  as  to  make  the  calculations  impractical, 
unless  an  effective  convergence  acceleration  scheme  can  be  found  (7:36). 

Accurate,  nonnegative  spatial  quadrature  schemes,  in  particular  the  exponential 
characteristic  (EC)  method,  are  effective  for  optically  thick  absorbing  regions  but 
with  unaccelerated  source  iteration  they  are  prohibitively  slow  to  converge  in 
thick  diffusive  regions. 

In  the  work  presented  here,  I  have  developed  new  algorithms  that  invert  the 
scattering  operator  in  each  cell  and  directly  solve  the  linear  system  of  coupled 
equations.  This  approach  eliminates  source  iteration  (SI),  per  se,  but  does 
require  iteration.  The  iteration  converges  cell  coupling  coefficients  which  depend 
on  the  angular  distribution  of  the  flux  and  (for  EC)  the  spatial  distribution  of  the 
source.  I  have  implemented  and  benchmarked  a  code  to  execute  one  dimension 
slab  geometry  particle  transport.  This  transport  is  efficient  in  thick  diffusive 
problems  for  the  two  non-negative  spatial  quadrature  schemes  tested.  The 
method  overcomes  the  inefficient  dependence  on  numerous  particle  flights  by  SI 
to  estimate  scattering  source  and  the  corresponding  lack  of  robustness  in  the 
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converged  SI  solution.  The  method  is  not  designed  as  an  accelerator  for  SI  but 
as  a  new  transport  method.  Since  the  new  method  converges  on  a  solution  by 
iteratively  seeking  improved  angular  and  spatial  distribution  information  I  refer 
to  it  as  the  angular  and  spatial  distribution  iteration  (ASDI)  method.  The  ASDI 
method  performs  comparably  with  the  SI  method  for  those  problems  that  don’t 
require  acceleration  and  rapidly  converges  in  problems  that  conventionally 
require  SI  acceleration.  This  is  particularly  important  in  the  thermal  neutron 
energy  range  from  about  0  to  1  eV  (^:83). 

In  systems  that  are  optically  thick  and  scattering  dominated,  particles 
undergo  many  collisions  before  being  captured  or  leaking  out  of  the  problem. 
Developing  a  practical  efficient  iterative  method  for  these  problems  is  of 
significant  practical  importance. 

Examination  of  the  one  group  time  independent  Boltzmann  Transport  Equation 
(BTE)  illustrates  the  problem.  A  one  group  particle  problem  in  planar  geometry 
can  be  expressed  as: 


A 


dx 


J  cr^  [x,ju'  ^  ju)i//(x,ju')  dju'  +Q^^\x,ju)  0  <  X  <  A, 


-1 


(1) 


r  (0,  a)  =  ( a)  +  Iq  ( A '  ^  a)  (0,  A ') 


//  >  0 


(2) 


r(^,A)  =  rr'*"^(A)  + j°j^^A'«i?(A'^A)r(^,A')  A<0  (3) 

where  x  is  the  position  coordinate;  /u  is  the  direction  cosine  of  the  angle  of  flight 
relative  to  the  positive  x-axis;  cr(.^)  is  the  total  cross  section;  cr^  (x,Q'-Q)  is  the 
scattering  cross  section;  g(x,//)  is  the  interior  emission  source,  ^(x,//)  is  the 
angular  flux  to  be  determined.  Equations  (2)  and  (3)  are  general  covering  all 
but  periodic  boundary  conditions  on  the  left  and  right  sides  with  appropriate 
choice  of  boundary  condition  i^ju'  ^  ju)  or  /r) . 
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By  defining  a  streaming  and  collision  operator  as 


I  =  (//|-  +  cr(x)),  (4) 

and  a  scattering  operator  as 

1 

5*  =  I  ,  (5) 

-1 

with  //^the  scattering  angle,  eqnation  (1)  is  written 

Ly/{x,ju)=  S\i/{x,ju)  +  Q‘^\x,ju).  (6) 

Generally  an  analytic  solntion  for  y/{x,iu)  is  not  possible.  Conventional  practice 
is  to  approximate  a  solntion  iteratively  nsing  Sonrce  Iteration  (SI).  The  SI 
scheme  is 


Ly/^‘-^^\x,fi)=  %g^^(x,//)  +  e£^(x,//) .  (7) 


Operationally  SI  works  as  ontlined  in  algorithm  I. 
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The  symmetric  relative  difference  function  £’§^0  is  introduced  in  Algorithm  1. 
This  function  has  all  the  properties  required  of  a  distance  function  (5:23),  as  is 
shown  in  appendix  A.  Thus,  the  real  numbers,  with  this  metric,  form  a  metric 
space.  This  metric  is  combined  with  the  vector  norm  |Hloo  to  measure  the 
maximum  distance  between  two  vectors,  as  is  explained  in  appendix  A. 

The  symmetric  relative  difference, y/^  ,  is  used  to 

determine  when  two  successive  flux  iterations,  (^x,ju)  and  y/^^^  [x,ju) ,  meet 

convergence  tolerance.  The  iteration  estimate  for  angular  flux  is  {x,ju)  ■  It  is 
the  angular  flux  due  to  particles  that  have  scattered  at  most  /  -1  times.  When 
particles  undergo  few  collisions,  the  SI  scheme  converges  rapidly.  However,  for 
problems  that  contain  diffusive  regions  that  are  optically  thick  and  scattering 
dominated,  SI  schemes  converge  slowly  and  may  converge  falsely.  If  scattering 
ratios  are  nearly  one  then  the  error  in  the  final  iterate  can  be  much  greater  than 
the  difference  (between  it  and  the  previous  iterate)  that  satisfied  the  preassigned 
convergence  criterion.  (/^:  10).  This  makes  it  difficult  to  determine  when  an 
iteration  scheme  is  suitably  converged  and  renders  an  SI  solution  unreliable. 

Slow  and  false  convergences  dictate  the  need  to  either  accelerate  SI  or  develop  a 
more  efficient  iterative  scheme. 
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Motivation 


Recently  Adams  and  Larsen  conducted  a  comprehensive  review  of  40  years 
of  methods  that  improve  iterative  transport  convergence.  This  work  outlines  four 
desirable  properties  for  fast  iterative  schemes  (4:139).  An  iterative  method 
should: 

1.  Converge  effectively  requiring  few  iterations.  Convergence 
effectiveness  is  typically  characterized  by  spectral  radius,  or 
equivalently  iteration  count.  Often  lower-  order  Sn  schemes,  or 
coarser  angular  refinement  are  proposed  to  accelerate  higher-order 
Sn  schemes.  Typically  these  lower  order  or  coarser  schemes  use 
fewer  unknowns  per  cell  and  require  fewer  transport  calculations,  or 
require  transport  calculations  that  are  computationally  cheaper 
than  the  higher-order  transport  scheme.  However  these  lower-order 
schemes  may  not  reduce  iteration  count  significantly  in  difficult 
problems. 

2.  Be  computationally  efficient.  Computational  efficiency  can  be 
measured  by  basic  memory  storage  requirements,  algebraic  cost  to 
implement  the  method  or  overall  computing  time.  If  the  equations 
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used  to  accelerate  a  transport  scheme  carry  the  same  number  of 


unknowns  or  use  computationally  expensive  cell  variables  the 
equations  may  yield  satisfactory  spectral  radius  yet  be  so  time 
consuming  to  solve  that  the  low  order  scheme  is  unacceptable. 

3.  Be  applicable  to  both  heterogeneous  and  homogeneous  problems. 
Recently  it  was  discovered  that  synthetic  acceleration  when 
extended  to  multiple  dimensions,  can  degrade  or  diverge. 

Transport  synthetic  acceleration  (TSA)  diverges  (5:15/16;  5:12/17) 
and  diffusion  synthetic  acceleration  (DSA)  degrades  from  a  spectral 
radius  of  1/3  in  a  homogeneous  material  to  a  spectral  radius  of  0.88 
(4:  139)  in  problems  with  periodic  material  interfaces.  Other 
acceleration  methods  are  expected  to  show  similar  degradation. 

The  transport  community  needs  a  scheme  that  maintains  a  small 
spectral  radius  and  is  cheap  computationally  in  both  slab  geometry 
and  multiple  dimensions. 

4.  Be  portable  to  parallel  systems  of  computers  and  not  degrade  in 

parallel  performance  as  the  number  of  processors  becomes  large. 

In  addition  to  the  four  properties  outlined  by  Adams  and  Larsen  I  add  a  fifth 
property.  An  iterative  method  should 
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5.  Be  robust.  By  this  I  mean  that  it  should  be  usably  accurate  for  the 


full  range  of  problems  that  we  want  to  solve.  This  accuracy  should 
not  be  achieved  by  operator  tuning,  fix-ups,  or  expert  system 
hybridization  but  by  aptness  of  the  algorithm  for  the  problem. 

This  property  will  be  referred  to  as  robustness. 

Synthetic  acceleration  and  quasidiffusion  techniques  have  been  applied  to 
particle  transport  in  problems  that  are  many  mean  free  paths  thick,  with 
scattering  ratios  close  to  unity  with  mixed  results.  In  slab  geometry  the  greatest 
improvement  in  efficiency  results  from  DSA  provided  that  the  diffusion  equation 
is  consistently  differenced  and  scattering  is  either  isotropic  or  weakly  anisotropic. 
I  sought  to  improve  the  iterative  convergence  rate  of  the  exponential 
characteristic  (EC)  method  as  developed  by  Mathews  and  Minor  in  1993(7) 
which  is  essentially  the  same  as  the  non-linear  characteristic  (NC)  method 
introduced  by  Waring,  Walters  and  Morel  in  1996  {8:  24-37).  Wareing  and 
Morel  subsequently  were  able  to  develop  an  effective  acceleration  method  in  slab 
geometry  for  their  NC  for  both  homogeneous  and  heterogeneous  materials  with 
scattering  ratios  of  one.  They  did  not  demonstrate  application  to  periodic 
material  interfaces,  and  they  did  not  report  compute  time  leaving  the 
computational  efficiency  of  their  acceleration  method  in  question.  Their  research 
has  not  been  extended  to  multiple  dimensions  (5:76).  Recent  research  has 
pointed  toward  the  failure  of  DSA  when  used  in  multiple  dimensions  on  problems 
with  periodic  interfaces.  Because  of  this  I  did  not  pursue  Wareing  and  Morel’s 
DSA  accelerator  for  use  with  EC.  Instead  I  sought  to  find  an  efficient  slab 
geometry  scheme  that  would  not  suffer  from  DSA-like  degradation  in  multiple 
dimensions.  My  motivation  for  doing  this  is  that  no  accelerated  source  iteration 
technique  to  date  is  unconditionally  stable,  yields  the  same  solution  as  the 
unaccelerated  source  iteration,  is  rapidly  convergent,  is  demonstrated  to  be 
computationally  efficient,  is  general  with  respect  to  geometry  and  can  be  applied 
to  the  EC  method.  The  research  effort  at  the  Air  Force  Institute  of  Technology 
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requires  efficient  solution  of  the  transport  equation  in  all  dimensions.  Mathews 
and  his  research  team  require  an  effective  (low  iteration  count),  computationally 
efficient  (fast  run  time),  robust  (reliably  accurate  without  user  intervention)  EC 
method  applicable  to  homogeneous  and  heterogeneous  materials  of  any 
configuration  or  geometry  in  order  use  the  EC  spatial  quadrature  on  a  wide  range 
of  transport  problems  of  interest  to  the  defense  community.  Further  it  is 
desirable  that  such  a  method  be  readily  adapted  for  use  with  parallel  computing 
environments. 


Goal  of  the  Research 


My  goal  was  to  develop,  implement,  and  evaluate  an  effective,  computationally 
efficient,  robust  method  for  solving  the  one  group,  slab  geometry  Boltzmann 
transport  equation  (BTE)  discretized  in  angle  and  space.  I  further  sought  to 
develop  a  method  that  was  general  with  respect  to  material  properties,  was 
readily  portable  to  parallel  computing  environments,  and  could  be  extended  to 
multiple  dimensions.  I  designed  the  method  primarily  for  use  with  EC  and  step 
characteristic  (SC)  spatial  quadratures  but  sought  a  method  that  could  also  be 
used  with  other  spatial  quadratures. 


Scope 
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I  derive  a  new  transport  method  that  explicitly  solves  for  infinite  particle  flights 
in  a  single  iteration.  The  method  conples  across  cells  in  space  nsing  a  two 
direction  angnlar  qnadratnre.  This  is  referred  to  as  the  global  space  solntion. 

The  method  farther  solves  particle  transport  within  each  spatial  cell  conpling  the 
N  directions  of  a  fine  angnlar  qnadratnre.  This  is  called  the  fine  angle  solntion. 
The  theory  of  both  the  global  space  solntion  and  fine  angle  solntion  is  introduced 
as  well  as  a  technique  to  combine  these  methods  in  an  iterative  scheme  that 
converges  flux  distribution  as  opposed  to  scattering  source  in  order  to  solve  the 
BTE  iteratively. 

The  method  is  implemented  and  tested  using  discrete  elements  one  group 
isotropic  average  cross  sections  with  EC  and  SC  spatial  quadratures.  The 
method  is  derived  in  a  way  that  generalizes  to  include  discrete  ordinates 
Legendre  moment  generated  cross  sections  (  ct/  ) ,  multigroup  anisotropic  cross 
sections,  and  other  positive  spatial  quadratures  such  as  linear  discontinuous  (LD) 
and  NC.  The  method  was  not  tested  with  spatial  quadratures  that  produce 
negative  fluxes  such  as  diamond  difference  (DD).  The  method  might  be 
expanded  to  include  these  spatial  quadratures  but  would  have  to  account  for 
negative  flux  values  in  calculation  of  flux  weights.  This  was  not  derived  or 
tested.  The  new  method  for  particle  transport  is  validated  by  comparison  with 
unaccelerated  conventional  SI  for  EC  and  SC.  The  symmetric  relative  difference, 
number  of  iterations  and  compute  time  for  the  two  methods  is  the  basis  of 
comparison. 

The  scope  of  the  test  problems  examined  is: 

•  Fixed  source,  sub-critical,  time  independent  systems 

•  Slab  geometry 

•  Single  group  isotropic  cross  sections 

•  Isotropic  emission  sources  uniform  in  each  cell 
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EC  and  SC  spatial  quadratures 


Although  the  method  was  examined  with  the  above  scope  it  was  derived  in  a 
general  way.  The  extension  to  multigroup  problems  is  immediate.  The  emission 
source  can  include  down  scatter,  up  scatter  and  fission  contributions.  Non¬ 
positive  spatial  quadratures  such  as  DD  will  require  new  algorithms  to  calculate 
flux  weights  from  negative  angular  fluxes. 

Assumptions  and  Limitations 


The  method  as  designed  and  tested  when  using  the  EC  spatial  quadrature 
and  discrete  elements  angular  quadrature  inherits  limitations  of  these 
approximations.  EC  is  not  currently  extendable  to  curvilinear  geometry. 

Discrete  elements  angular  quadratures  lead  to  numerical  approximations  that  are 
not  the  same  as  the  diffusion  approximation. 


Approach 


An  explicit  system  of  equations  for  the  angular  and  spatially  discretized  BTE 
using  a  general  linear  spatial  quadrature  is  derived  in  N  directions.  This  system, 
although  illustrative,  is  impractical  to  solve.  A  similar  system  using  only  two 
directions  is  introduced.  A  closed  form  solution  of  this  two  direction  fully 
spatially  coupled  system  is  derived  which  is  correct  for  transport  with  only  two 
directions.  The  concepts  of  transport  coefficients  and  flux  weights  are 
introduced.  Transport  coefficients  define  the  relationships  of  incoming  angular 
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fluxes  and  source  emissions  to  outgoing  angular  fluxes  within  a  spatial  cell.  Flux 
weights  facilitate  the  projection  of  the  system  of  equations  between  multiple 
directions  and  two  directions.  Then  they  are  used  to  calculate  a  closed  form 
solution  using  a  coarse  angle  approximation  which  effectively  solves  for  flux 
across  the  entire  problem  space.  The  cell  edge  flux  calculated  in  this  way  is 
correct  if  the  flux  weights  used  to  collapse  the  transport  coefficients  are  correct. 
The  coupling  of  the  exact  N  direction  transport  method  (using  approximate  edge 
flux  values)  with  the  equivalent  2  direction  collapsed  transport  coefficient  method 
(also  using  approximate  edge  flux  values)  is  then  introduced.  An  iteration 
scheme  that  produces  progressively  more  accurate  flux  weights  which  are  used  to 
compute  progressively  more  accurate  transport  coefficients  is  then  discussed. 

This  method  effectively  amounts  to  iteration  on  transport  coefficients.  The 
method  differs  significantly  in  concept  and  execution  from  accelerated  SI.  For 
instance,  accelerated  source  iteration  follows  the  logic  of  algorithm  2. 


Intialize  Qw 


Calculate  Ql  y/ 


Apply  Boundary  Conditions 


SolveLy/^^^^^  (^x,ju)  =  Q^^^  [x,ju)  for 


Apply  accelerator  correction 

TPjn  rl  rl  n  c _ { i/r  ,A  1/r 


r‘nici-\y/yv*rF/:>irir‘/y  /  nl/yv*nirir‘/y 


The  method  proposed  in  this  research  follows  the  logic  of  Algorithm  3 


Initialize  cell  coupling  coefficients  with  coarse  angle  approximation 
Use  2  direction  SC  cell  coupling  coefficients  to  generate 

Do 

Solve  for  fine  angle  flux  within  cells 
Generate  improved  low-resolution  cell  coupling  coefficients 
Solve  for  fine  angle  flux  within  cells 
Use  low  resolution  edge  flux  to  improve  high  resolution  edge  flux 

I  r  .,\  !  r\ 

Algorithm  3  is  implemented  in  Fortran  95  and  benchmarked  against 
nnaccelerated  sonrce  iteration  which  is  also  implemented  in  Fortran  95. 

The  ASDI  method  is  then  extended  to  the  EC  spatial  qnadratnre  by  adding  an 
iteration  loop  to  find  the  sonrce  distribntion  parameter,  fl ,  bnt  otherwise  retains 
the  same  logic  as  algorithm  2.  It  is  also  written  in  Fortran  95  and  is 
benchmarked  against  an  nnaccelerated  sonrce  iteration. 
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II.  Solution  of  the  Discrete  Ordinates  {S„)  Transport  Equations 


Solution  of  the  monoenergetic  BTE  introduced  in  Chapter  1  requires  treatment  of 
the  space  and  angular  variables  through  a  number  of  discretization  techniques  that  yield 
simultaneous  equations.  The  discrete  ordinates  method  is  a  widely  used  method  for 
obtaining  numerical  solutions  to  the  integrodifferential  form  of  the  transport  equation. 

It  is  the  method  used  in  this  research  to  discretize  the  Boltzmann  transport  equation  in 
space  and  angle.  A  brief  discussion  of  space  and  angle  discretization  follows. 

In  one  group,  slab  geometry  the  scattered  source  on  the  right  hand  side  of  equation  (1) 
may  be  replaced  by  a  Legendre  expansion  (2:13,36,117): 

00 

Y,(2h+l)P{p)cT,^^  0<x<X, 

k=0 

where  are  the  Legendre  moments  of  the  scattering  cross  section  and  ^4  the 
Legendre  moments  of  the  angular  flux: 

1 

<Pk{^)=  \  (9) 

-1 

The  discrete  ordinates  approximation  consists  of  requiring  equation  (8)  to  hold  only  for 
a  number  (N)  of  discrete  directions  then  applying  a  compatible  quadrature 
approximation  to  the  flux  moments  in  each  of  these  directions.  The  weighted 

quadrature  used  for  the  flux  moments  has  the  form 

N 

A{x)='Z^n'Pk{Mn')w{x,Mn')-  (10) 

n'=\ 

Substituting  the  weighted  quadrature,  at  N  distinct  angles,  with  a  Legendre  expansion 
truncated  at  a  finite  number  of  polynomials  (K)  in  equation  (8)  results  in 
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(11) 


N  K 

Z  Z  (2^  +  l)^^i  {^)Pk  {^n  )Pk  (/“«')  0  <  X  <  X 

«'=1  k=G 

The  coupling  of  the  system  of  equations  occurs  on  the  right  hand  side.  The  left  hand 
side  represents  a  set  of  angularly  independent  first  order  differential  equations.  Defining 
cr„  to  be  the  ordinate  to  ordinate  scattering  cross  section: 

^n,n 

(^)  =  Z  (2^  +  (^)^k  {Mn’)Pk  i^n  )  ,  (12) 

k=0 

facilitates  writing  equation  (11)  as 

^n^^»^^Ka{x)if/„ix)=  ^  a,^,^^{x)w„,i//„,{x)  +  Qf  {x)  .  (13) 

n’=l 

Equation  (13)  represents  a  system  of  N  coupled  differential  equations. 

The  method  I  have  developed  is  dependent  on  representing  the  discrete  ordinates 

equations  in  the  form  of  equation  (13)  with  the  coupling  of  the  system  taking  place  on 

N 

the  right  hand  side  through  the  term  ^  <j^  ,  (x) .  How  the  cross  section 

n'=\ 

data  weights  and  ordinate  direction  cosines  are  arrived  at  is  not  important  to  my  work. 
Other  methods,  such  as  that  introduced  in  2003  by  Gerts  and  Mathews,  which  calculate 
cross  sections  using  piecewise  averages  and  discrete  elements  {10),  may  be  used  if  the 
discrete  direction  flux  equations  can  be  represented  in  a  form  analogous  to  equation  (13) 
A  brief  discussion  of  their  work  follows. 

Dropping  spatial  dependence,  the  scattered  source  term  on  right  hand  side  of 
equation  (1)  is  written 

Q,=  j  (14) 

4;t 

where  Q  is  a  unit  vector: 
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Q  =  {Sin {^9) Cos Sin {^9) Sin {^(o), Cos {^9)) , 


(15) 


and 

ju  =  Cos{^9^  (16) 

The  discrete  elements  method  approximates  scattered  sonrce  within  an  element,  ,  by 

a,,=  |a(n)dn.  (17) 

An„ 

and  angnlar  flnx  within  an  element  by 


¥n'=  j  (18) 

In  the  above  eqnation  AQ„  is  the  Cartesian  prodnct  of  the  angnlar  interval  from  [0,2;r) 
and  the  //  interval  containing  the  ordinate  .  The  nnion  of  the  p.  intervals  covers  the 
range  [—1,1]  with  no  overlap.  An  element  of  solid  angle  on  the  snrface  of  the  nnit 
sphere  in  ID  slab  geometry  is: 

2n 

j  dCl=  j  did^dco.  (19) 

AQ„ '  A//„ '  0 

Eqnation  (18)  is: 

2n 

=  j  dju  ^  dco  \j/[ju,(i>)  ^  (20) 

ISjUn'  0 

Eqnation  (17)  is: 

Ik 

Ss„  =  j  j  (21) 

A/7„  0 
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and  equation  (14)  is 


In  In 

Qsn  -  j  dju^dco^  j  J//' j  (22) 

0  n'  0 

The  order  of  summation  and  integration  in  equation  (22)  may  be  rearranged  resulting 
in: 

2n  2n 

a,=I  I  dju^dco  I  J//' I  .  (23) 

n'  ls.n^  0  A//„'  0 

In  slab  geometry  y/i^iu\o}'^  has  no  co'  dependence: 

=  (24) 

Substituting  equation  (24)  into  equation  (20)  results  in: 


¥n'=  j 

A//„. 

Substituting  equation  (25)  into  equation  (23)  results  in: 


(25) 


2n  2n 

«'  A//„  A//„. 


(26) 


0  0 


Within  each  angular  element,  the  flux  is  approximated  as  isotropic: 


Wn' 


(27) 


Substituting  equation  (27)  into  equation  (26)  and  rearranging  the  order  of  integration 
results  in 
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The  integral  of  equation  (28)  defines  the  element  to  element  scattering  cross  section 

T  I  2^  j  '^TC 

=  I  ‘i/1  I  do>'iT^({/i',o>'}.{/i,a}).  (29) 

1..  aI  i  i 

Using  this  discrete  elements  angnlar  qnadratnre, 

Qsn=lL^^n'^n¥n'  ■  (30) 

n' 

The  element-to-element  scattering  cross  section  is  analogons  to  the 

ordinate-to-ordinate  scattering  cross  section  in  eqnation  (13).  With  discrete  elements, 
there  are  no  qnadratnre  weights  per  se.  These  are  replaced  by  the  sizes  of  the  elements 
Aju„ ,  and  are  treated  implicitly  throngh  the  calcnlation  of  .  Bothe  angnlar 

qnadratnres  have  the  same  form  if  weights  w„  =1  are  introdnced  in  eqnation  (30) .  Gerts 
and  Mathews’  work  has  the  advantage  of  representing  cross  sections  with  non- negative 
valnes  in  sharp  contrast  to  certain  cross  sections  arrived  at  throngh  Legendre  expansion. 
My  derivation  and  implementation  is  done  with  discrete  elements  cross  sections  because 
they  are  non-negative.  I  will  use  the  term  ordinate  in  my  work  since  it  is  more  common 
in  the  transport  community  but  I  intend  for  the  term  to  refer  to  both  ordinate  in  the 
discrete  ordinates  sense  and  element  in  the  discrete  elements  sense.  The  numerical 
results  presented  in  chapter  3  and  chapter  5  used  discrete  elements. 

Because  the  coupling  of  the  system  represented  by  (13)  occurs  on  the  right  hand 
side,  it  is  useful  to  combine  the  right  hand  side  as  a  single  term  (x) : 

=  +  (31) 

where 
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(32) 


N 


S„(X)=  X 

n'  =  \ 


The  source  of  scattered  particles  is  (x)  and  the  source  of  emitted  particles  is  (x) . 
Equation  (13)  is  then 


Mn 


^Wnix)  ^  ^ 
dx 


{x)\l/n{x)  =  g„(x) 


n  =  \,...N. 


(33) 


This  results  in  N  coupled  differential  equations  for  N  ordinates. 

The  system  represented  by  equation  (33)  is  easily  solved  once  spatially 
discretized  using  an  iterative  scheme  like: 


+  =  df\x) 


n  =  \,...N 


(34) 


where  is  an  iteration  index.  This  is  called  source  iteration.  If  one  starts  with  a  guess 
of  zero  for  scalar  flux  the  first  iteration  yields  the  uncollided  flux,  the  second  iteration 
yields  the  uncollided  plus  first  collided  flux,  the  next  iteration  yields  the  uncollided  plus 
first  and  second  collided  flux  and  so  on.  In  order  to  calculate  a  numerical  solution  for 
equation  (34),  it  is  first  spatially  discretized.  The  indexing  for  spatial  discretization  is 
shown  in  Figure  I. 
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X  =  0 


x  =  X 


Xq 


Xi 

2 


Xi 


Xi  Xs 

2  2 


Xj-i 

X;_|  Xj 


^I+\ 


xj-; 


1 

2 


Figure  1  Spatial  Index 

In  the  spatial  discretization  shown  in  Fignre  1,  cell  edge  qnantities  carry  half¬ 
integer  indices  and  cell-average  qnantities  carry  integer  indices.  The  total  nnmber  of 
mesh  cells  is  denoted  by  /.  On  the  left  and  right  sides  are  phantom  cells,  denoted  with 
0  and  I-l-l  respectively.  These  phantom  cells  are  a  convenient  acconnting  constrnct. 
They  produce  neither  intrinsic  nor  scattered  source.  They  have  no  thickness  so  they 
produce  no  particle  losses.  They  do  however  contain  edge  flux  information  at  the 
problem  left  and  right  boundaries. 

The  first  step  in  the  spatial  discretization  of  equation  (34)  is  to  spatially 
integrate  it  over  a  cell  to  obtain  a  balance  equation  for  the  spatial  cell.  In  the  i**"  cell, 
as  indicated  by  Figure  1,  the  left  edge  is  at  x._i  while  the  right  edge  is  at  x.  i  .  The 

I  2  '  +  2 

thickness  of  cell  i  is  Ax^-  =x.  i  -x._i  .  We  presume  that  material  discontinuities  are  also 

i  +  2  '“2 

cell  edges  so  that  cr(x)  is  a  constant  cr^- ,  within  each  cell.  Angular  flux  on  the  right  and 

left  cell  edges  are  y/  1,1//  1  respectively.  The  subscript  n  denotes  the  ordinate  and 

n,i+—  n.i— 

2  2 

the  subscript  i  denotes  the  mesh  cell.  Integration  of  the  of  the  flux  in  the  cell 
results  in 
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(35) 


-^(r  .  \  -¥  .  \  )+  <^i¥„.i  =  Q„, , 

At-  n,i+—  n.i—  ’  n,i 

AAA;  .  2  ’2 


Where  Q  .  is  the  average  source  in  cell  i: 


n^i 


X.  1 

“ 


e,„.=J,‘7*2.w- 


(36) 


The  average  source  in  cell  i  consists  of  average  scatters  and  average  emissions 


Q  ■  ~ 

n^i  n^i  n^i 


(37) 


Again  the  average  scattered  source  in  cell  i,  5*^  ,  ,  is  found  by  integrating  the  scattered 


n^i 


source,  5'„(x)  ,  over  the  cell: 


n^i 


(38) 


The  average  intrinsic  source  in  cell  , ,  is  also  found  by  integrating  the  intrinsic 
emission  source,  (x) ,  over  the  cell: 


X.  I 

f  i& 

Jx  I  ^ 

'A 


(39) 


In  practice  the  average  scattered  source  is  not  calculated  from  equation  (38)  but 
from  scattering  cross  sections  and  average  angular  flux.  The  average  angular  flux  in  cell 
i  is: 


(40) 

A 

The  scattering  cross  sections  from  ordinate  n'  to  ordinate  n  are  assumed  to  be  constant 
for  cell  i  (as  was  done  for  total  cross  section)  and  are  denoted  as  cr„  ,  .  Hence,  the 

'  •^n,nj  ' 

average  scattered  source  of  equation  (38)  can  be  written 


2-8 


(41) 


N 


Using  equation  (41)  and  equation  (37)  to  substitute  for  the  right  hand  side  of  equation 
(35)  yields  a  system  of  equations  that  can  be  solved  analytically  if  another  equation  is 
introduced  and  one  of  the  fluxes  (incoming,  outgoing,  or  average)  is  known.  If  the 
incoming  flux  is  known  for  each  direction  in  a  cell,  the  average  and  outgoing  fluxes  are 
the  system  unknowns.  Two  equations  are  needed  for  each  direction  but  only  the 
balance  equation  in  each  direction  has  yet  been  introduced.  To  complete  the  system  of 
equations  an  additional  equation  for  each  direction  is  needed.  This  additional  equation 
is  known  as  the  auxiliary  equation.  Given  the  balance  equation  and  an  auxiliary 
equation  for  a  cell  with  an  N  direction  angular  refinement,  it  is  straightforward  to  solve 
for  the  exiting  flux.  This  can  be  done  for  every  cell  in  a  spatial  approximation. 

The  conventional  scheme  sweeps  through  the  cells  in  the  direction  of  motion  of 
the  particles.  Consider  an  example  with  quadrature  points  indexed  in  order  of 
decreasing  direction  cosine  but  not  necessarily  symmetric  about  ju  =  0  .  Let  ju^  through 
ju^^he  positive  (rightward)  and  through  ju^  be  negative  (leftward).  The 

solution  process  for  a  multicell  spatial  discretization  begins  by  solving  for  exiting  edge 


flux  y/^  +i/_i  average  flux  ,  where  /  denotes  the  total  number  of  mesh  cells 

and  Njf  + 1  denotes  a  leftward  direction  whose  direction  cosine  is  smallest  in  magnitude 
among  the  direction  cosines  of  the  quadrature  set.  The  incoming  edge  flux  at  the  right 


face  is  y/^  +i/+i  ■  known  from  the  boundary  condition  at  the  right  edge.  With  a 

vacuum  or  source  boundary  these  values  are  explicitly  known  but  with  a  reflective 
condition  the  incoming  edge  flux  values  are  set  to  the  appropriate  outgoing  edge  flux 
values  using  the  latest  iteration  estimate  for  these  fluxes.  After  solving  for  outgoing  flux 


at  edge  y/^  +i/_i  value  is  used  as  the  input  for  the  adjacent  cell  then  y/^  +i/_3 

calculated.  The  solution  process  commonly  referred  to  as  a  sweep 
proceeds  across  each  cell  until  the  left  boundary  is  reached.  The  process  is  repeated  for 
all  remaining  fluxes  with ju^<0 .  When  the  left  face  is  reached  for  each  ordinate  the 


This  value  is  used  as  the  incoming 
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edge  flux  for  the  far  left  cell.  Outgoing  edge  flux  ( y/^  3 )  and  average  angular  flux 
are  calculated.  Outgoing  edge  flux  from  cell  1  ( y/  3 )  is  used  as  input  for  the  second 

^  R  ’“2 

cell.  This  procedure  continues  through  the  mesh  successively  solving  for  the  average 
and  outgoing  fluxes  in  each  cell  until  the  right  boundary  is  reached.  The  process  is  then 
repeated  for  all  remaining  ordinates.  An  iteration  is  complete  when  fluxes  for  all 
ordinates  are  calculated  in  all  cells.  Upon  completion  of  an  iteration,  the  source  term  is 
re-calculated  from  equation(37)  and  the  sweep  process  begins  again. 

An  alternative  to  source  iteration  is  to  consider  the  entire  spatial  and  angular 
mesh  as  a  system  of  equations  with  outgoing  cell  edge  fluxes  and  cell  average  fluxes  as 
the  unknown  variables.  These  variables  are  calculated  by  multiplying  cell  incoming 
fluxes,  average  intrinsic  sources,  and  average  scattered  sources  by  coefficients  that  are 
determined  by  the  spatial  quadrature  (spatial  differencing)  scheme.  These  spatial 
quadrature  coefficients  are  scalar  quantities  that  are  calculated  from  cell  cross  section 
<J^ ,  particle  direction  ju^  ,  and  cell  thickness  .  Certain  spatial  differencing  schemes, 
such  as  EC,  also  use  first  moments  to  calculate  them.  I  will  discuss  EC  in  chapter  4. 

The  discrete  ordinates  system  of  equations  represented  by  equation  (33)  can  be 
written  as  the  contribution  of  cell  incoming  edge  flux,  scattered  source  and  emission 
source  multiplied  by  a  transport  coefficient  to  cell  outgoing  edge  flux  is 


^  .1  S  \y/ 


n,i  -^n,!  /  n,i- 


+  ’  Mn  ^  ^  Sa^  , ,  ^ 


+  Koe„  ,  ( ,  Av,. ,  Sa^  .  ,  .  1 . 


>0, 


(42) 


n,i —  \  /  n,i-\ — 


+  ,■  (  ’  Mn  ^  ^  Sa^  , ,  ,  A  , 


+  Koe„  ,  (  ,  Mn  ^  ^  ^A„  , .  , 


<0. 


(43) 


The  spatial  quadrature  also  provides  coefficients  for  determining  the  cell  average  fluxes: 
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¥n,i  =  ^^4  ^  ^  ,  S  W 


n,i— 


+  Kas„  ,■  (  ’  Mn  ^  ^i  ^  Sa„  , ,  . 


+  Kae„  ,  ’  Mn  ’  ’  Sa„  , ,  >^x„  ,• 


/«„  >0, 


y^n,i  =KaI„, 

’  2 

+  Kas,j  (ct.  ,  /<„ ,  Ai.  ,  . ,  S,__ , )  Sj^. 

+  ^AE„ ,  (^i ,  M„ .  AVi ,  Sj. , ,  S,. , ) 


/!„  <0. 


(44) 


(45) 


For  compactness  of  notation  the  transport  coefficients  will  be  written  without 
their  arguments,  unless  these  arguments  are  necessary  for  discussion  or  clarity.  The 
spatial  quadrature  coefficients  will  be  written  as  Kqj^,,  Kq^^,,  , 

Kai^  . ,  Kas„  i  ’  I  ■  These  spatial  quadrature  coefficients  can  be  obtained  from  any 
spatial  discretization  (spatial  quadrature).  They  are  a  convenient  representation  of  the 
coupled  differential  equations  of  the  discrete  ordinates  equations.  These  coefficients  are 
derived  for  an  SC  spatial  quadrature  in  the  next  section.  SC  is  chosen  in  this  research 
because  it  is  a  positive  method,  like  EC,  but  is  linear  which  simplifies  implementation. 
Successful  SC  implementation  does  not  guarantee  extension  to  EC  but  was  used  to  test 
method  implementation  before  applying  it  to  the  more  sophisticated  EC  spatial 
quadrature  which  uses  first  moments. 


SC  Transport  Coefficients 


In  this  section,  formulas  for  the  step  characteristic  (SC)  spatial  quadrature 
coefficients  are  derived.  SC  approximates  the  scattering  and  emission  sources  as  being 
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uniformly  distributed  in  a  cell.  The  flux  in  the  cell  is  obtained  as  characteristic 


solutions  along  each  direction  of  the  angular  quadrature.  This  is  evaluated  at  the 


outflow  face  and  is  averaged  over  the  cell  to  obtain  these  coefficients. 


Dividing  equation  (33)  by  results  in: 


dx 


(46) 


Let  the  optical  thickness  from  a  point  x  to  a  point  x  along  the  direction  of  the 
ordinate  be  (^x,x^: 


X  <  x'  <  X  , 


(47) 


If  £7  is  a  constant  between  x  and  x ,  then 


(48) 


T  1  X  X  I 

Now  introducing  the  integrating  factor  e  ”  ^  '  into 


/  Ar^\  1  1  •  •  1  dur  (x)  , 

equation  (46)  and  replacing  — ^ -  with  - - results  m 


5x 


dx 


dx 


(49) 


Noting  that 


A. 

dx 


(X)/"  (*■*))=  /"  t'*)  (,v)/”  (*■*) , 

)  dx 


(50) 
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and  integrating  equation  (49)  from  x  to  x  produces 


,  ,  T„(x,x)  r~\  T„(x,x)  rx  dx'  ^  r„(x’,x) 

il/„{x)e"^  ’ -i//Jx)e'^^  ^  =  L - Q„{x')e"^ 

^  ^  Mn 


(51) 


Dividing  equation  (51)  by  the  integrating  factor  and  rearranging  terms  results  in 


¥ni^)  =  ¥[x)^ 


jjni.x,x)  r„(x,x))  j^rxd^  q  ^r„(x’,x)-r„(x,x) 

Mn  ^  X  1 1 


(62) 


T  \  X  x\ 

where  e  ’  'is  1.  If  //„  >  0  the  flux  is  entering  from  the  left  cell  edge  and  x  =  x^.  If 
//„  <  0  the  flux  is  entering  from  the  right  cell  edge  and  x  =  x^  .  This  leads  to 


y/n{x)=y/{xL)^ 


T^n{x,XL)  J_  ^  Q  ^rni.x',Xi)-T„{x,Xi) 

^^n  JxL 


//„>0,  (53) 


wni=‘)=¥(xRi  -  yj"'— e„(x') 

”  f^n  •^x  u 

H-n 


Hn  <0.(54) 


In  the  i“"  cell  as  indicated  by  Figure  1  x^  corresponds  to  x._i  while  x^  corresponds  to 

'  2 

X,  Further  if  Ax  =  x.  j^-x.  j^  is  sufficiently  small,  cr(x)  can  be  assumed  to  be  a 
1+2  1  +  2  '  2 

constant  cr^- ,  the  average  over  the  cell.  Equation  (48)  is  then 


(55) 


This  is  readily  expressed  in  terms  of  the  optical  thickness  of  cell  i  in  direction  n  ,  which 
I  denote  as 


Mn 


(66) 


Thus, 
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x,x,  , 

V  '27 


=  \S, 


r  A 

x-x.  I 
V  ‘~2J 


I  Ajc,. 


(57) 


and  in  a  similar  way 


x,x,  I 

V  '47 


=  \S. 


X  I  -X 

'4 


n,i\  Ax 


(58) 


Snbstitnting  equation  (57),  equation  (58),  and  the  source  integral  as  given  by  equation 
(31)  these  cell  flux  equations  are: 


f  ) 

\^n,i 

X-X  1 

|l  'aJ 

AX; 

(x-x') 

Ax,- 


(x-x') 

Ax; 


(59) 


r  A 

X.  I  -X 
V  2  y 

Ax/ 


. 


|(x'-x) 


(x'-x) 

Ax,- 


(60) 


/'»<0 


Physically  these  equations  determine  the  flux  y/^iXi)  in  a  cell  as  the  superposition  of  the 
flux  entering  at  a  cell  edge  and  the  flux  due  to  source  production  from  scatters  and 
emissions.  In  the  above  equations  S„(x')  is  the  source  of  scattered  particle  emissions 
and  £'„(x')  is  the  source  of  intrinsic  particle  emissions  at  position  x'  into  ordinate  n.  A 
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characteristic  scheme  introduces  approximation  by  assuming  a  convenient  form  for  the 
scattered  source  distribution  then  computes  the  exact  angular  flux  corresponding  to  this 
assumed  source  distribution.  A  step  characteristic  spatial  quadrature  assumes  that  the 
source  distribution  is  constant  in  cell  i  in  ordinate  .  This  constant  source 
distribution  averaged  over  a  spatial  cell  is 


n,i 


rVi  dx 

•'b-i  Av 
2 


S„ix). 


(61) 


In  practice  is  not  calculated  from  equation  (61)  but  is  obtained  from  available  cell 

information  using 


N 


n’=l 


n,i 


¥n\v 


(62) 


Setting  x  =  x.  i  in  equation  (59)  produces  the  flux  exiting  the  right  edge  of  spatial  cell  z 
given  by 


.  1  =1^  .  1  e 

n,i+—  n.i— 

2  2 


r 

■Sa  ^ 


i+j  dx' 


kl 


r 

+Ea  A  ‘ 
1 


z+^  dx' 


1 

f  ) 

e 

X.  1  -x’ 

Ik  J 

Axi 

f  ] 

1 

e 

X.  1  -x' 

Ik  J 

Axj 

4  Wn\ 


>0. 


(63) 


Similarly  setting  x  =  x.  i  in  equations  and  (60)  produces  flux  exiting  the  left  edge  of 
spatial  cell  i  given  by 


Equation  (60)  solved  for  flux  exiting  the  left  edge  of  spatial  cell  i  is: 
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W  .  \=¥  .  le'  ’ 
n,i —  n,i-\ — 

2  2 


r  \ 


X.  1 

-Sa  . f 


Wn\ 


f  ’ 

+Ea  .[  ‘ 


x’-x.  1 

1 

ISXj 

f  ) 

x’-x  1 

l^«,/ 1' 

V  ij 

Ax; 

(64) 


4  \^n\ 


Mn<^- 


X  -X,  1 

Introducing  the  following  change  of  variables  to  m  '  =  — -r — -  in  equations  (63)  and  (64) 


results  in 


.  1  =1^  .  1 
nj+—  nj-— 

2  2 

\^n,i 

e  '  ' 

o  ^  “k 

+  Sa  .  1 - re  1 

\m,\ 

"’’\Cdu' 

Jo 

,  Z7  Ax  -  £ 

+Ea  .  — re  ' 

w  / 

’  \Mn\ 

n,i  — 

2 

nA-\ — 

2 

+  *^^4  -1 

kf" 

'^E A  . 

(65) 


>0 


(66) 


<0. 


Evaluation  of  the  integrals  in  equations  (65)  and  (66)  yields  a  positive  scalar  quantity 
for  each  transport  coefficient.  Practical  evaluation  of  these  integrals  requires  using 
exponential  moment  functions  as  introduced  by  Mathews  and  Minor  (ii:169).  The 
exponential  moment  function  of  order  zero  with  one  argument  is: 

1 

31^  (x)  =^— — .  (67) 

0 
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The  exponential  moment  function  of  order  n  with  one  argnment  is: 


{-xt)  _ 


Using  eqnation  (67)  to  evalnate  the  integrals  of  eqnation  (65)  results  in 


¥  .  i=¥  . 

n,i-\ —  n,i  — 

2  2 


Ax  - 


Ax  -\e„  i 


Hn  >0. 


Noting  that 


-X 


Mq{x), 


results  in  rewriting  equation  (69)  as 


.  1  =1^  ,  1  e  '  ’  ' 

nj-\ —  ltd  — 

1  2 

H-n 


i^(K,|) 


Mn  >0- 


Using  equation  (67)  to  evaluate  the  integrals  of  equation  (66)  results  in 


¥  .  i=¥  .  le 
rid —  nd-\ — 

2  2 


\j^n 

Ax 


r^n 


Mn  <0- 
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Equations  (71)  and  (72)  are  in  the  form  of  equations  (42)  and  (43).  The  spatial 
quadrature  coefficient  for  the  outgoing  flux  due  to  the  incoming  flux,  K^j ,  is  obtained 
from  yhe  first  term  on  the  right  of  equations  (71)  and  (72): 


K, 


oi 


n,i 


(73) 


The  second  transport  term  on  the  right  hand  side  relates  the  contribution  of  flux  within 
cell  scattering  to  the  flux  exiting  an  edge: 


Mn 


(74) 


The  third  term  on  the  right  hand  side  relates  the  contribution  of  cell  intrinsic  emissions 
to  the  flux  exiting  an  edge: 

=  — ^(K.,|)-  (76) 

Mn  ^  ^ 


Note  that  for  SC,  =KQg^, . 

The  coefficients  for  the  cell  average  fluxes  are  developed  next.  Equations  (59) 
and  (60)  averaged  over  a  spatial  cell  are: 


.=r^: 

Jx.  I 


^■'-"4  dx 


—  ¥  .  1  e 

Av  n,i-- 

2 


x-x  1 

I  '2 
'  Ax 


X  I  r  I  l(^“^') 

^  r  ;+4  dx  rx  dx'  — 

+  ^  ^ 


I 


f 

C,!  Jx.  j 


^  dx 


i+-k  ax  rx  dx' 

I  _ z? 


l(x-x') 
'  Ax 


i  Ax-’Vi|//„| 


(76) 


>0, 
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f  1+^  dx 
M=l  .  1 


X  1  -X 

I 

'  Ax 


+  —I 

^tlA  '  ‘  ' 


r- 

'2 


i+^  dx  r  *+]■  dx'  -'• 


li„ 


tA  fVi  dx'  -  ^n 


|(x'-x) 
I  Ax 


|(x'-x) 
I  Ax 


(77) 


x-x.  1  X  -X.  1 

Using  the  change  of  variables  u  =  — ,  and  w '  =  ■  ^ 


Ax 


//„  <0. 


for  eqnation  (76)  then 


rearranging  terms  resnlts  in 

fl 

'n,i  =\J^V^ 

’  JO  n,, 

pi 

r  du 

Jo 


+  Sj  .  — A  due  I  ’  I  Jm'  e'  ’ 

JO 


+  E. 


Ax  fl 


(78) 


"’‘kro  -"o 


//„  >0. 


Using  eqnation  (67)  to  evalnate  the  first  term  of  eqnation  (78)  and  evalnating  the  right 
most  integral  of  the  second  and  third  terms  of  eqnation  (78)  results  in 


¥n,i  (K,l|) 


+  S 


Ax 


+  E 


\/An\  r»j 
Ax 


1  r  ^ 

du 

■A  Jo 


A„i 


kl  hi\io 


1  f  ^ 

du 

Jo 


1  -  e 


e  1^”’' 

hi 

y 

e~hj\ 

y 

(79) 


/!„  >  0. 


Using  equation  (67)  to  evaluate  the  remaining  integral  in  the  second  and  third  terms 
results  in 
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¥n,i  ._1-H  (K,y|) 


+  s. 


Ax 


+  E, 


Wn\ 

Ax 


i-^(h,/|)' 


^n,i 


i-^(k/|) 


Mn  >0- 


(80) 


Substituting  equation  (68)  where  appropriate  in  equation  (80)  results  in 

¥n,i  ._1-H  (K,^|) 

n,i  ^ 

’  \Mn\  ^  ^ 

x-x.  1  x’-x.  1 

^7  ^"*’7 

Using  the  change  of  variables  u  - — ,  and  u'  =  — for  equation  (77)  then 
rearranging  terms  results  in 


Using  equation  (67)  to  evaluate  the  first  term  of  equation  (82)  and  evaluating  the  inner 
integrals  of  equation  (82)  results  in 


’  1 

r  ^ 

du 

Jo 


Ax  I 


g  r«,;r_g 


(83) 


+  E, 


Ax  1 

r  ^  1 
du 

f  \  \  1  1  A 

g  r  g  r”>6 

En 

Jo  1 

V  y 

En  <0- 
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Using  equation  (67)  and  equation  (70)  to  evaluate  the  remaining  integral  of  equation 
(83)  results  in 


¥n,i  = 


Ax 

f  1-^(1- 

'n,i  J 

1 

^fi,i 

y 

Ax 

'n,i 

1) 

1 

) 

Substituting  equation  (68)  where  appropriate  in  equation  (80)  results  in 

n,i+—  ''I 

+  5. 


\Mn 

Ax 


Mn  <0- 


(84) 


(85) 


Equation  (81)  and  equation  (85)  are  again  in  the  form  of  equation  (44)  and 
equation  (45).  The  fourth  transport  coefficient  relates  the  contribution  of  cell  entering 
edge  flux  to  cell  average  flux  it  is 

(86) 

The  fifth  transport  coefficient  relates  the  contribution  of  cell  scattered  flux  to  cell 
average  flux.  It  is 


K 


AS, 


(87) 


The  sixth  transport  coefficient  relates  the  contribution  of  cell  intrinsic  emissions  to  cell 
average  flux.  It  is 
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(88) 


=1^31,  (K, I)  ■ 

\M'n  \ 

Note  that  for  SC  =Kap 

The  spatial  quadrature  coefficients  in  this  section  are  derived  for  SC.  A  similar 
procedure  is  followed  with  other  spatial  quadratures  which  result  in  different  formulas 
for  the  spatial  quadrature  coefficients  that  also  can  be  put  in  the  form  of  equations  (42) 
and  (43).  The  use  of  these  coefficients  in  the  ASDI  method  is  general  and  applicable  to 
coefficients  derived  from  any  spatial  quadrature. 


Explicit  Solution  of  Transport  Equations  for  N  Directions  Coupled  in  Space 


The  cell  transport  equations;  (63),  (64), (76),  (77),  written  in  vector  notation  are 


(89) 

1^4-  ini  ^ A  '^^AEi  ^A  ’ 

(90) 

SAi  =s^.  y/^, , 

(91) 

where  bold  type  denotes  a  matrix,  denotes  a  vector,  and  i  denotes  a  cell  number.  If 
the  total  number  of  ordinates  is  denoted  by  N  ,  the  ordinate  set  is 
I  //j ,  V  •  •  5  Mnp  ’  ’  ■  ■  ■  Atv  }  •  Denoting  the  rightward ,  >  0  ,  ordinates  as 
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yjUi,jU2,--- ^  Nji  is  the  number  of  rightward  directions.  The  leftward,  ju„<0, 
ordinates  are  .  The  number  of  leftward  directions  is 

Nl=N-Nj,.  (92) 


The  inward  flux  vector  in  cell  i  ( )  has  N  components 


V^in.-  = 


^Nn+U+\ 


(93) 


The  top  set  of  angular  fluxes  in  equation  (93)  represent  the  flux  in  at  the  left  edge  of 


cell  i  streaming  in  the  rightward  direction, ._i  |,  and  the  bottom  set 
represent  flux  in  at  the  right  streaming  in  the  leftward  direction  in  cell  i 


The  outward  flux  vector  in  cell  i  ( W out- )  ^  components 


^  out; 


(94) 
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The  top  set  of  angular  fluxes  in  equation  (94)  represent  the  flux  out  at  the  right  edge  of 


cell  i  streaming  in  the  rightward  direction,|p^^  | ,  and  the  bottom 

represent  flux  out  at  the  left  edge  of  cell  i  streaming  in  the  leftward  direction 


set 


Equation  (93)  can  be  written  even  more  compactly  as 


y^in.  = 


y' in  n. 


y^i, 


(95) 


where 


y^inR 


NR,i 


,■  1 


(96) 


and 


^4,  = 


¥ 


NR+\,i+\ 


(97) 


Similarly  equation  (94)  can  be  written 


¥  out: 


¥^ 


out  R. 


¥, 


out  R. 


(98) 


where 
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w 


OUtj^. 


•  1 


W«,!  + 


(99) 


and 


W 


OUti- 


(100) 


The  average  intrinsic  emission  vector  in  cell  i  {Ea^  )  has  N  components.  Written  in 
compact  notation  it  is 


where 


and 


^^R,i 

^^L,i 


EAu,  = 


(101) 


(102) 


(103) 


The  coefficient  matrix  denoting  the  contribntion  to  the  ontgoing  flnxes  from  the 
incoming  flnxes  in  cell  i  is  a  diagonal  matrix  of  size  NxN: 
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0 


0 


0 


Ko/.  - 


K, 


Oh 

0 

0 

0 


1,; 


K, 


oi 

0 

0 


2,i 


0 

0 


0  K, 


OI 


N,i 


(104) 


The  other  transport  coefficient  matrices  ’^ae^  ’  have  the  same 

form  as  equation(104).  They  are  diagonal  matrices  of  their  corresponding  spatial 
qnadratnre  coefficients. 

The  matrix  contains  the  element-to-element  scattering  cross  sections,  or  the 
ordinate-to-ordinate  scattering  cross  sections  and  weights: 


i  ■ 

^^N^2,i 

(106) 

^S2^N,i'^2  ■ 

The  following  matrices  are  nsefnl  for  stripping  the  rightward  or  leftward  components  of 
the  vectors  introdnced 


Ir  - 


0  0 


0  0 


0  I 


Nl 


(106) 


(107) 


where  and  I^^are  identity  matrices  of  dimension  Nj^  and  Ni  respectively,  and  0 
represents  a  matrix  of  zeros  of  appropriate  dimensions.  As  an  example  of  how  these 
matrices  are  nsed  I  will  rewrite  eqnation  (98)  in  terms  of  the  ontgoing  flnx  from  the  two 
cells  adjacent  to  cell  i.  Using  eqnation  (106)  and  equation  (107)  results  in 
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(108) 


W  OUtji^ 

Inr  0 

W  outjf. 

+ 

1 - 

o 

> 

o  < 

1 _ _ _ 

W  OUtjl. 

W  OUti^ 

0  0 

W  OUti. 

VJ 

1 - 

o 

1 _ 

An  explicit  solution  of  the  discrete  ordinates  equations,  (89)  and,  (90)  is  possible.  First 
consider  cells  that  are  not  on  either  the  left  or  right  exterior  boundaries  so  that  particles 
stream  into  the  cell  from  the  adjacent  cell  on  each  side  and  particles  stream  out  of  the 
cell  into  those  cells.  Substituting  equation  (91)  into  equation  (89)  and  (90)  results  in 


Wouti  =Ko/,.  ’  (109) 

=^AIi  ¥ini  ^Si  +^AEi  •  (110) 


The  formal  solution  of  equation  (110)  for  y/ results  in 


=(/-K 


ASi 


)-iK 


Ali  --^ASi  -e. 


(Ill) 


Substituting  equation  (111)  into  equation  (109)  with  some  rearrangement  results  in 


^ouU-^Oh  +^05  S  (/-K^_5  S  )  W 


I  \ 


i  I  '  ‘"i 


+  |Ko£,  +Ko5.S^,(/-K^^,  Ea,. 


(112) 


This  may  be  written  more  compactly  as 


V^outi  =^OIi  Win.  +^OEi  EAi  , 


(113) 


where 


^oii  -  \  ^oii  +  ^oSi  I  ’ 


(114) 


and 
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I  will  refer  to  the  elements  of  the  matrices  ^oii  ™OEi  transport  coefficients. 

The  flux  entering  at  the  left  of  cell  i  moving  rightward  is  ( )  .  The  flux  exiting 

V  '"i  )n'<Nji 

at  the  right  of  cell  i  moving  rightward  is  ( )  ■  The  fraction  of  flux  entering  at 

the  left  ( //  >  0 )  and  leaving  at  right  ( //  >  0  )  is  ( m^r.  ) 

V  '  fn,n' 

Using  equation  (108),  equation  (113)  can  be  written 

^outi  =^OIi  {^R)¥outi_^  +^OIi  +^OEi  ^ •  (116) 

Equation  (113)  expresses  the  outgoing  edge  flux  for  a  mesh  cell  )  in  terms  of  cell 

emission  (  Eai  )  and  outgoing  edge  flux  from  the  two  adjacent  cells.  The  coefficient 
matrix  vaqj.  determines  the  contribution  to  cell  outgoing  flux  from  cell  incoming  flux. 
The  coefficient  matrix  ^qEi  determines  the  contribution  to  cell  outgoing  flux  from  cell 
intrinsic  emission. 

Next  consider  cells  at  the  left  or  right  boundaries.  If  a  cell  is  at  an  exterior 
boundary,  particles  stream  into  the  cell  from  one  adjacent  cell  and  particles  behave  in 
accordance  with  a  boundary  condition  on  the  other  side.  The  reflection  condition  for 
this  boundary  condition  for  any  discrete  ordinates  set  can  be  represented  as  a  matrix. 
This  matrix  is  denoted  by  cv,^  on  the  left  boundary  and  cv,^  on  the  right  boundary. 

Each  element  of  the  matrix  represents  the  fraction  of  flux  from  an  outgoing  ordinate 
that  is  reflected  into  an  incoming  ordinate.  Defining  these  boundary  matrices  as 
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1 

o 

...  0 

ar 

^1,Nr+1 

ar 

^\,N 

•  o 

...  0 

aj 

^Nr,Nr+\ 

ar 

^Nr,Nr 

■  •  o 

...  0 

■  ■  o 

■  ■  o 

1 

o 

...  0 

•  o 

1 

..  o 

0 

0 

0  ••• 

0" 

0 

0 

0  ••• 

0 

Oi,^  - 

^^Nr+1,1 

^Nr+\,Nr 

0  ••• 

0 

_  ^Rn,i  ■■■ 

^^N,Nr 

0  ••• 

0 

(117) 


(118) 


facilitates  specification  of  all  but  periodic  boundaries.  In  block  form  these  boundary 
matrices  are: 


where 


(119) 


Otr  = 


aj 


a. 


4,Afs+l 


ai 


aj 


^N,N 


(120) 


and 


OLr  - 


0 

Q.^ 


where 


(121) 
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(122) 


Oi,R  - 


a 


w,i 


Consider  a  problem  without  incident  flux,  cell  1  which  is  at  the  left  boundary  has  as  an 
incident  flux  on  the  left  side  the  flux  resulting  from  the  above  boundary  condition 
rather  than  the  inward  flux  denoted  in  equation  (93).  Using  the  above  boundary 
condition  this  inward  flux  at  the  left  boundary  ( y/-^^ )  is  written: 


N 


n'=NR  +  \ 


,  1 
’2 


¥in, 


N 


S 

n'=NR  +  \ 


CCt 


^n\NR 


.  i 

’2 


w 


Nr+\, 


3 

2 


(123) 


Equation  (123)  can  be  rearranged  as 


y'in,!  = 


N 


n'=N R+\ 


N 


y  ocj V ,\ 

^  _  ^n’,NR  ^  n',^ 


n'=N R+l 


0 

0 


+ 


0 


0 


(124) 


Equation  (124)  represents  the  inward  flux  of  the  left  boundary  as  the  contribution  to 
the  cell  flux  from  the  left  side  boundary  reflection  and  particles  streaming  in  from  the 
right  adjacent  cell.  In  compact  matrix  notation  this  is 
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(125) 


V'in,  • 

In  any  spatial  discretization  scheme  adjacent  cells  are  conpled  at  the  edge  that  is  shared 
by  the  two  cells.  This  adjacent  cell  conpling  for  cells  one  and  two  for  the  leftward 
streaming  flnx  is  expressed  as: 

^L¥in,=^L¥out^-  (126) 


With  this  expression  (125)  can  be  written 

¥in,=^L¥out,+^L¥out2-  (127) 

Snbstitnting  equations  (127)  into  equation(113),  and  applying  the  left  side  boundary 
condition  results  in 


¥outy  =  “o/i  ¥out,  +  “o/i  Iz  Wout2  +  (128) 

Rearranging  terms  in  equation  (128)  results  in 

(l  -  mo/j  ^4  •  (129) 

The  right  most  ( )  cell  in  the  spatial  domain  is  handled  in  the  same  way  resulting  in: 

(l  -  =  moij ^RWoutj_i  +  ^4  •  (130) 

The  system  represented  by  equation  (116)  has  2N  equations  in  2 4/ unknowns 
corresponding  with  interior  cells.  Equations  (129)  and  (130)  each  represent  N 
equations  and  N  unknowns  corresponding  to  the  two  cells  on  the  boundaries.  There 
appears  to  be  more  unknowns  than  equations  until  one  considers  the  connectivity  of  the 
mesh  cells.  This  connectivity  results  from  the  edges  shared  by  adjacent  cells.  For  an 
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interior  cell  the  incoming  edge  flux  is  the  outgoing  edge  flux  from  the  adjacent  cell. 
This  leads  to  N  connectivity  equations  for  the  (T2)  interior  cells 


ini  outi_i 

/  =  {2,. 

(131) 

L-Wim  =^L¥outi^^ 

z={l,...,7-l}. 

(132) 

Note  that  the  incoming  edge  fluxes  on  any  cell  edge  can  be  written  as  the  sum  of  the 
outgoing  edge  fluxes  from  two  adjacent  cells. 

On  the  left  boundary  in  cell  1 ,  flux  enters  from  edge  z  =  1  /  2  and  from  edge 
z  =  3/2 .  The  incoming  flux  in  cell  1  at  edge  z  =  1/2  is  the  outgoing  flux  at  this  same 
edge  reflected  back  in  the  incoming  direction.  The  incoming  flux  in  cell  1  at  edge  i  = 
3/2  is  the  outgoing  edge  flux  from  cell  2.  In  cell  1  there  are  N/2  connectivity  equations 
and  N/2  boundary  equations.  Likewise  the  right  edge  incoming  flux  is  the  outgoing  flux 
reflected  at  the  boundary  at  edge  z  =  7  +  1/2  and  from  the  adjacent  upstream  cell  i  =  I- 
1/2.  In  cell  I  there  are  N/2  connectivity  equations  and  N/2  boundary  equations.  The 
boundary  cells  contribute  N  connectivity  equations  and  N  boundary  equations  to  the 
spatially  coupled  system.  Collectively  there  are  2NI  unknowns  and  2NI  equations. 

A  particle  transport  system  of  equations  that  connects  all  the  interior  cells  with 
the  boundary  cells  and  couples  all  directions  represents  fully  coupled  discrete  ordinates 
transport.  This  angular  and  spatially  coupled  system  of  transport  equations  can  be 
written  in  the  following  way.  Let  '¥in  represent  incoming  flux  vectors  for  the  cells  in  a 
spatial  mesh  of  I  cells  with  a  phantom  cell  on  the  left  side  and  a  phantom  cell  on  the 
right  side.  This  vector  has  2  N  (I+I)  components  and  is  written 


Let  out  represent  outgoing  flux  vectors  for  the  same  mesh: 


^out  = 


outQ 

outi 


y^  out  I 
y^  outj_^Y 


(134) 


Let  Ea  represent  average  intrinsic  emission  vectors  for  the  same  mesh,  it  is  written: 


Ea 


Eaj 

0 


(135) 


The  zeros  at  the  top  and  bottom  of  the  emission  vector  of  equation  (135)  represent  left 
and  right  phantom  cells  which  do  not  emit.  In  equation  (133)  ,  W outi^  represent  the 

flux  entering  or  exiting  a  phantom  cell  at  the  left  boundary,  Wmi^^  ■>  y'outj^i  represent  the 
flux  entering,  or  exiting  a  phantom  cell  at  the  right  boundary.  Flux  enters  and  leaves 
these  phantom  cells  through  a  single  edge.  Using  equations  (93)  through  (100)  these 
arrays  can  be  written 


y^itiQ 


y^L. 


'OUt\ 


(136) 


V 


OUtQ 


y^R, 


-mi 


(137) 
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(138) 


'«/+! 


y^R. 


outj 

0 


0 


y^L. 


‘OUtj-^\ 


(139) 


where  is  a  vector  of  zeros  of  length  or  .  Because  the  incoming  edge  flux  for 
interior  mesh  cells  is  the  outgoing  edge  flux  from  an  adjacent  cell,  can  be  expressed 
as  a  combination  Using  the  matrices  specified  by  equation  (106)  and 

equation  (107)  is  written 


li? 

01 

y^in,  = 

0 

- 1 

1 - 

1 

o 

+ 

1 _ 

The  phantom  cell  incoming  flux  at  the  left  edge  is  written 


1 — 

o 

o 

1 _ _ 

1 

o 

1 _ 

y^iiiQ 

- 1 

z 

o 

_ _ _ 1 

y^T 

^  out\ 

(140) 


(141) 


The  phantom  cell  incoming  flux  at  the  right  edge  is  written 


and  is  written 


In«  0 

- 1 

1 

>2 

O 

_ 1 

0  0 

1 

O 

_ 1 

li? 

0  “ 

- 1 

1 

>3 

+ 

_ 1 

0 

1 - 

1 

T 

i _ 

The  phantom  cell  outgoing  flux  at  the  left  edge  is  written 


(142) 


(143) 
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(144) 


outQ 


In«  0 

0  0 

1 

o 

_ 1 

The  phantom  cell  outgoing  flux  at  the  right  edge  is  written 


0 

0  ■ 

0 

^  OUtjj^Y 

0 

1 - 

1 

a' 

1 _ 

(145) 


Using  equations  (140)  through  (145)  the  global  spatial  vector  Tm  can  be  written  in 
terms  of  the  global  spatial  vector  Toti/  : 


0 

¥p 
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^out\ 
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0 
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0 

0 

0 

0 

0 

^  Ln,jU 

^in\ 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

ouiY 

^Rinl 

0 

0 

ItVr 

0 

0 

0 

0 

0 

0 

0 
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0 

¥p 

^out2 

0 

0 

0 

0 

0 

0 
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^Nl 

0 

0 

0 

0 

¥r 

^out2 

^Rinj_^ 

0 

0 

0 

0 

Inr 

0 

0 
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0 
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0 
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0 
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^Nr 

0 

0 

0 

0 

0 

¥r 

^Rinj 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

0 

¥p 

^outj 

0 

0 

0 

0 

0 

0 

0 

0 

ItVr 

0 

0 

0 

¥i 

^OUtj 
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0 

0 

0 
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0 

0 

0 

0 

0 

0 

0 

0 

0 

y^T 

(146) 


Two  of  the  equations  leading  to  the  system  denoted  by  equation  (146)  carry  no 
information.  These  correspond  to  the  phantom  cell  edge  flux  corresponding  to  the  edges 
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that  are  not  shared  with  cell  1  or  cell  1.  Noting  that  the  first  and  last  rows  and  the 
second  and  next  to  last  colnmn  are  zero  in  the  matrix  of  equation  (146)  we  can  remove 
these  equations  from  the  system  since  they  carry  no  information.  Doing  this  results  in 
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Defining  P^^^as  the  matrix  of  equation  (147) 


.(147) 
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0  0  0  0  0  0  •  ••  0  0  0  0 
0  0  0  0  0  0  •  ••  0  0  0  0 

0  0  0  0  0  0  •••  0  0  0  0 

0  0  0  0  0  0  •  ••  0  0  0  0 

0  0  0  0  0  0  Inl  •  ••  0  0  0  0 

P/o=  ,  (148) 

0  0  0  0  •  ••  Inj^  0  0  0  0  0  0 

0  0  0  0  •  ••  0  0  0  0  0  0 
0  0  0  0  •  ••  0  0  Inj^  0  0  0  0 

0  0  0  0  •  ••  0  0  0  0  0  0  1^^ 

0  0  0  0  •  ••  0  0  0  0  0  0 


equation  (147)  can  be  written  more  compactly  as 

^in  =  P/o  ^out.  (149) 

Note  that  is  a  permutation  matrix  that  reorders  W out  to  become  'Fm  .  As  such,  its 
inverse  is  its  transpose. 

Global  spatially  coupled  coefficient  matrices,  denoted  as  and  ,  are 

assembled  from  the  cell  coefficient  matrices  of  equations  (114)  and  (115).  These 
matrices  are  ordered  so  that  their  first  and  last  rows  correspond  with  phantom  cells  and 
their  other  rows  correspond  with  interior  cells 
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m„  m„, 

,Nr  ,1  ,Ni  ,1 

°  ,Nj{  ,1  <^^Nl  ,Ni  ,1 

m„,  m„, 

o^Nji,Nji,2  o^Nji,Ni,2 
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UlnP  UlnP 

0 

(151) 

Using  equations  (114),  (135),  (150),  and  (151),  the  spatially  coupled  system  becomes: 

'^out  =Mqj  'Yin  +Mqj-  Ea  .  (152) 
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Equation  (152)  represents  the  transport  equations  explicitly.  It  accounts  for  all  particle 
penetration  and  scattering.  Using  equation  (149)  and  combining  'Fom/ coefficients  results 
in 


(l-Moi  ^io)^out  =Moe  Ea  .  (153) 

The  matrix  pjis  reordered  to  result  in  a  banded  diagonal  matrix  of  minimum 

bandwidth.  The  reordering  matrix  which  achieves  this  minimum  bandwidth  structure  is 
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Applying  P^to  equation  (153)  and  solving  for  out  yields  the  explicit  solution  of  the 
fully  coupled,  in  both  angle  and  space,  system  of  transport  equations 

=  (Pa  (l  -  Mo/  P,o ))"'  P/,  Mof  £./  .  (155) 

Edge  flux  computed  from  equation  (155)  yields  an  explicit  solution  for  the  discrete 
ordinates  equations  without  iteration.  It  accounts  for  particle  transport  within  a  cell 
and  particle  transport  across  cells.  Unfortunately,  the  2 A/ system  which  results 
from  A  directions  and  /  cells  is  impractical  to  invert  or  solve,  particularly  if  the  spatial 
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and  angular  discretization  are  fine.  This  is  why  source  iteration  (SI)  does  not  attempt 
to  solve  this  system.  Instead  SI  approximates  the  solution  by  iteration  on  the  scattering 
source. 


Angular  and  Spatial  Distribution  Iteration 


I  did  not  attempt  to  solve  the  system  of  equations  represented  by  equation  (I55)in  this 
research.  Rather,  I  investigated  an  efficient  iteration  scheme  that  was  different  from 
source  iteration.  I  solved  equation  (113)  with  a  refined  angular  quadrature  locally 
within  each  individual  spatial  cell  using  an  estimate  for  incoming  edge  flux  and  its  cell 
transport  coefficients.  I  then  estimated  particle  transport  among  cells  across  the 
problem  space  using  equation  (155)  with  a  coarse  angular  quadrature  and  transport 
coefficients  that  are  correct  for  the  approximate  edge  flux  used.  I  developed  and 
examined  a  method  that  couples  local  within  cell  scattering  using  a  refined  angular 
quadrature  with  across  cell  particle  transport  using  transport  coefficients  obtained  from 
a  coarse  angular  quadrature  in  a  unique  way  that  retains  angular  resolution.  The 
concept  is  similar  to  2  direction  synthetic  acceleration  but  differs  from  it  in  two 
important  ways: 

1.  The  proposed  method  does  not  use  transport  sweeps.  The  proposed 
method  does  not  approximate  the  error  in  the  flux  iterate  to  correct  a 
transport  sweep. 

2.  The  method  estimates  the  solution  directly  with  a  coarse  angular 
quadrature  and  a  fine  angular  quadrature  linked  with  transport 
coefficients. 
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The  method  uses  the  coarse  angular  quadrature  to  directly  solve  equation  (155).  This 
provides  edge  fluxes  (Tm  )  that  are  approximate.  These  edge  flux  values  are  used  in 
equation  (113)  to  gain  a  better  estimate  of  transport  coefficients  {Kq^ ^  .,K^g ^  .,Kqj ^  . , 
Kjj ),  in  a  way  not  yet  discussed,  in  order  to  obtain  a  better  estimate  of 
T  in  for  a  two  direction  spatially  coupled  angular  quadrature.  In  my  method 
4 /equations  are  solved  in  the  spatially  coupled  coarse  angle  routine  or  4A^  equations 
are  solved  in  the  fine  angle  within  cell  routine  as  compared  with  the  2N I  equations  of 
the  fully  coupled  system,.  By  de-coupling  space  and  direction  the  system  of  equations  is 
reduced  to  either  four  times  the  number  of  ordinates  in  the  angular  quadrature  or  four 
times  the  spatial  cells.  This  is  far  less  than  the  number  of  ordinates  times  the  number 
of  spatial  cells  required  for  the  fully  coupled  system. 

Unique  aspects  of  this  research  are  the  use  of  the  flux  distribution  along  cell 
edges  to  yield  progressively  better  edge  flux  transport  coefficients  and  of  iteration  on 
angular  and  spatial  flux  distributions.  The  method  effectively  amounts  to  improvement 
of  transport  coefficients  and  iteration  on  angle  and  spatial  flux  distributions.  The 
method  eliminates  the  need  for  modeling  neutron  flights  by  transport  sweeps  as  is  done 
with  SI.  The  method  significantly  reduces  transport  iterations,  even  when  scattering 
ratios  are  nearly  1.  The  method  produces  reliably  accurate  answers  even  when  neutron 
flights  are  nearly  infinite.  Discussion  of  the  fine  angle  and  coarse  angle  solution 
methods  and  their  coupling  follows. 


Explicit  Solution  of  a  Two  Direction  Transport  Problem  Coupled  in  Space 


Although  it  was  not  practical  to  solve  the  discrete  ordinates  system  of  equations 
with  many  directions  coupled  across  space,  it  is  practical  if  there  are  only  two  directions 
to  couple  across  space.  Consider  a  slab  geometry  transport  problem  with  two 
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directions.  Following  the  vector  notation  of  the  previous  section  the  cell  transport 
equation  is  given  by 

=Moe.£.4  ,  (156) 

where  Mqj  is 
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and  Mq£  is 
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Pis  given  by  equation  (148)  and  by  equation  (154)  with 
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The  components  of  the  spatially  coupled  coefficient  matrix  ( )  are  given  by 
equations  (114)  and  (115).  The  diagonal  coefficient  matrices  used  in  these  equations  are 


K 


oii 


K, 


Oh,i 

0  Ko,^. 


(161) 


K 


OSi 


K, 


OSi,i 

^  ^OS2^i 


(162) 


^Ai..  = 


K 


^h,i 

0  K 


AI 


2,i 


(163) 


K 


ASi 


K 


0 


(164) 


The  scattering  cross  section  matrix  is  a  dense  matrix 


(165) 


If  there  are  /  cells  in  a  spatial  partition  the  arrays  of  incoming  fluxes,  outgoing  fluxes 
and  intrinsic  emissions  are 
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The  transport  coefficient  entries  of  the  left  hand  side  matrix  of  equation  (167) 

( ^oii  1  ,■  ’  ^oii  2  /  ’  ^oix  2  ,■  ’  ^oi2  2  i )  describe  the  fraction  of  entering  flux  in  a  direction  (left 
or  right)  that  leaves  the  cell,  having  scattered  any  number  of  times  within  the  cell.  A 
particle  enters  a  cell  from  one  of  two  possible  directions.  Having  entered  a  cell  in  an 
ordinate,  a  particle  then  either  transmits  to  the  other  side  in  that  ordinate  or  is 
reflected  backward  into  the  opposite  ordinate.  Because  of  this  it  is  useful  to  think  of 
the  transport  coefficient  entries  as  reflection  or  transmission  coefficients.  The 
transmission  coefficient  describes  the  fraction  of  particles  that  enter  a  cell,  scatter  any 
number  of  times,  then  exit  that  cell  in  the  original  direction  of  its  flow  while  a  reflection 
coefficient  describes  the  fraction  of  particles  that  enter  a  cell  scatter  any  number  of 
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times  then  exit  that  cell  in  the  opposite  direction  to  its  original  flow.  Defining  reflection 


coefficients  as: 


’  (168) 

h,i  =  ^Ol2,2i  ’  (^®9) 


and  transmission  coefficients  as 


"^0/1,2.  . 

(170) 

moi  2^1  ’ 

(171) 

eqnation  (167)  is  written 
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If  the  spatial  quadrature  conserves  particles  and  c  <  1 ,  then  0  <  ;■  <  1 ,  0  <  ;■  <  1 ,  and 

;■  <  1 .  If  c  =  1  then  +  =  1 .  If  the  boundary  conditions  of  the  system  are 

vacuum,  grey  or  albedo,  then  the  system  is  diagonally  dominant  when  c  <  1  and  weakly 
diagonally  dominant  if  c  =  1 .  Because  of  this,  any  ID  transport  problem  with  a  unique 
non-trivial  steady  state  solution  results  in  at  least  a  weakly  diagonally  dominant  system 
with  positive  values  on  the  diagonal  and  negative  values  off  the  diagonal.  Therefore  the 
ID  slab  geometry  system  with  two  ordinates  is  readily  solved.  If  an  angular  quadrature 
with  two  ordinates  adequately  described  the  angular  flux  distribution  of  a  problem  we 
would  be  able  to  calculate  by  inverting  the  matrix  of  equation  (172)  and  solve  the 
problem  directly. 

Unfortunately,  a  two  ordinate  angular  approximation  is  most  often  inadequate. 
For  example,  particles  that  enter  a  mesh  cell  parallel  to  the  +x-axis  are  more  likely  to 
exit  a  mesh  cell  in  the  +x  direction  then  particles  that  enter  a  mesh  cell  nearly 
perpendicular  to  the  +x  axis.  The  5*2  angular  quadrature  does  a  poor  job  of  modeling 
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directional  dependence  for  this  type  of  particle  transport.  Because  of  this,  an  5*2 
angular  quadrature  is  often  an  inadequate  approximation  for  angular  flux.  I  sought  to 
overcome  the  weakness  of  this  coarse  approximation  in  angular  distribution  yet  take 
advantage  of  its  ability  to  fully  couple  cells  in  space  without  numerous  particle  flights.  I 
speculated  that  because  the  coarse  angular  quadrature  estimated  flux  correctly  in 
magnitude  these  flux  values  could  be  used  as  reasonable  estimates  for  cell  entering  flux. 
This  required  development  of  a  projection  scheme  between  coarse  and  fine  angular 
quadratures.  I  also  speculated  that  the  approximate  edge  flux  values  obtained  from  the 
spatially  coupled  coarse  solution  could  be  angularly  coupled  to  the  flux  locally  within  a 
spatial  cell  in  order  to  approximate  angular  flux  distribution  on  cell  edges.  Further,  I 
speculated  that  the  approximate  angular  flux  distribution  obtained  this  way  could  be 
used  to  improve  the  two  ordinate  transport  coefficients.  I  expected  the  improved  two 
ordinate  coefficients  to  in  turn  improve  the  approximation  for  cell  incoming  flux.  I 
sought  to  use  the  computational  efficiency  of  the  coarse  angle  spatially  coupled  system 
and  the  fine  angle  local  space  systems  with  an  iteration  scheme  that  effectively  coupled 
the  two  schemes  through  transport  coefficients. 


Solution  of  an  TV  direction  Transport  Problem  with  Approximate  Edge  Flux. 


Explicit  solution  of  equations  (113)  through  (115)  is  possible,  even  for  very  fine  angular 
quadratures,  if  entering  flux  )  is  known.  Equation  (113)  may  be  used  to  solve  for 
exiting  flux  for  any  resolution  angular  quadrature  without  creating  a  system  that  is 
impractical  to  solve  since  it  is  done  for  only  a  single  cell.  The  inverse  matrix  of 
equations  (114)  and  (115),  models  scattering  within  a  cell  effectively 
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because  it  effectively  accounts  for  an  infinite  number  of  neutron  flights.  This  can  be 
seen  by  expanding  the  matrix  in  a  geometric  series: 
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(173) 


The  term  j  in  this  geometric  series  produces  the  flux  of  neutrons  that  have 

scattered  k  times  within  the  cell  (without  leaving  the  cell).  The  first  term  ,1,  accounts 
for  the  flux  of  neutrons  that  have  not  scattered  since  entering  the  cell.  Thus,  the  series 
accounts  for  all  numbers  of  scattering  events  within  the  cell  for  each  particle  entering 
the  cell  or  emitted  within  the  cell.  If  the  incoming  flux  is  correct  the  outgoing  flux  can 
be  explicitly  calculated.  In  application  y/-^  is  not  a  known  quantity.  It  is  an  estimate. 
This  estimate  can  be  improved  by  using  equation  (113).  Given  estimated  incoming  fine 
angle  flux,  equation  (113)  is  used  to  calculate  fine  angle  outgoing  flux.  This  is  done  for 
each  spatial  cell.  After  each  cell  outgoing  flux  is  calculated  boundary  conditions  are 
applied.  Outgoing  fluxes  calculated  in  this  way  are  incoming  fluxes  for  adjacent  cells 
with  fine  angle  resolution.  This  leads  to  an  improved  incoming  flux  estimate.  I 
examined  whether  edge  flux  produced  with  a  two-ordinate  angular  quadrature  from 
equation  (172)  could  be  used  as  input  for  a  fine  angle  flux  solution  produced  by 
equation  (113)  which  in  turn  could  produce  adjusted  transport  coefficients.  I  then 
examined  whether  the  adjusted  transport  coefficients  then  improved  the  approximate 
flux  of  the  coarse  angle  estimate.  This  provided  an  iteration  scheme.  I  further 
examined  the  convergence  properties  of  this  scheme.  The  combination  of  the  fine  angle, 
and  coarse  angle  methods  is  discussed  next. 
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Combining  High  Angular  Resolution  Within- Cell  Transport  with  Low  Angular 
Resolution  Spatially-Coupled  transport 


The  previous  sections  described  a  method  for  obtaining  exact  within  cell  solntions  for 
edge  flnx.  Continning  the  notation  developed  in  these  sections,  is  given  by  eqnation 
(93),  is  given  by  eqnation  (94),  E^,  is  given  by  eqnation  (101).  The  scattering 
cross  section  matrix  is  given  by  eqnation  (105).  The  coefficient  matrices 

,K^^,  ,K^^,  ,^ae^  diagonal  matrices  of  the  form  given  by  equation 

(104).  The  array  of  cell  average  flux  is  written 


y^Ni^+\,i 


(174) 


¥N,i 


With  this  notation  the  transport  system  is 

outi  V^itij  +^05;  ^5*.  V^Aj  ^Ai ,  (175) 

y^Aj  ^^Ali  V^inj  +^^5,  ^5.  y^ Aj  '^^AEi  '  (176) 

The  flux  of  particles  is  transported  through  the  cells  of  a  spatial  mesh  as  though 
the  particles  were  moving  in  representative  directions.  In  discrete  ordinates  each 
component  of  the  flux  vector  describes  the  flow  of  particles,  through  a  surface  patch  of 
solid  angle  on  the  unit  sphere.  Figure  2  shows  the  flux  exiting  a  cell  edge  in  the  left 
( //  <  0 )  and  right  ( //  >  0 )  hemispheres  for  both  coarse  and  fine  angular  refinements.  A 
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2  Directions 


Cell  Edge 


^2  =  400 


4  Directions 


^3  =  100 


^4  =  300 


coarse  angnlar  qnadratnre  of  two  directions  is  depicted  in  the  top  set  of  hemispheres 
and  a  finer  angnlar  qnadratnre  of  fonr  directions  is  depicted  in  the  bottom  set  of 

Figure  2:  Angular  Flux  Exiting  a  Cell  Edge  through  2  Ordinates  and  4  Ordinates 


hemispheres.  The  notation  \j/  denotes  a  flnx  associated  with  a  coarse  angnlar 
qnadratnre.  The  symbol  y/  denotes  a  flnx  from  a  finer  qnadratnre.  Elements  one  and 
fonr  are  the  caps  of  the  hemispheres  and  elements  two  and  three  are  the  bands  closes  to 
the  cell  edge  denoting  that  these  flnxes  are  in  fact  inclined  at  a  steeper  angle  to  the  x 
axis  which  passes  throngh  the  tip  of  the  polar  cap. 


The  snbscript  m  denotes  an  ordinate  in  a  two  direction  angnlar  qnadratnre  whose 
snrface  elements  are  the  two  hemispheres  in  Fignre  2.  The  total  flnx  exiting  a 
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hemisphere,  y/^  ,  is  calculated  by  summing  all  the  finer  angular  quadrature  fluxes 


exiting  the  same  hemisphere. 
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In  a  similar  way  the  average  flux  in  a  spatial  cell  in  a  hemisphere  denoted  as,  y/^  ■ ,  is 


and 
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y^2,i=  Z  ^n,i  /  =  {1,---,I}.  (180) 


Further  the  average  particles  scattered  per  volume  in  the  solid  angle  whose  average 
particle  direction  is  in  one  of  the  hemispheres  of  Figure  2  is  denoted  ■  It  is 

calculated  by  summing  the  particles  per  solid  angle  exiting  the  associated  hemisphere 
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The  average  particles  per  volume  emitted  in  the  solid  angle  whose  average  particle 
direction  is  in  one  of  the  hemispheres  of  Figure  2  is  denoted  ■  It  is  calculated  by 

summing  the  particles  per  solid  angle  exiting  the  associated  hemisphere 

£iu=i;£^„,  (183) 

n=\ 

<  =  (184) 

n=Nji 

As  in  the  previous  sections  denotes  the  number  of  rightward  flux  ordinates 
andA^  the  number  leftward 


N  =  Nj^+Ni^. 


(185). 


Defining  ■,¥ out-  fhe  coarse  angular  quadrature  inbound,  outbound,  and  cell 

average  angular  flux  arrays  in  cell  i,  they  are  written  as 
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Defining  the  array 


1  0 
0  1 


(190) 


where  the  rows  of  I2  have  N  elements  consisting  of  or  ones  or  zeros.  This 
matrix  facilitates  writing  coarse  qnantities  as  the  snm  of  fine  qnantities.  These  coarse 
flnxes  in  compact  form  are 


Win-  -  h,N  Win-  ’ 
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W outi  ^2, A'  W outi  1 

(192) 

W A^  =  I2,7V  W A^  ■> 

(193) 

EAi  =  l2,7V  EAj  . 

(194) 

The  transport  equations  for  these  two  directions  are  given  by  equations  (109)  and  (110). 
Operating  on  these  equations  with  I2  results  in: 


y^outi  -  Win-  '^\N^OSi'^Si 

WAi  =I2,7v(k^/.  J  +  I2,7V  ^5',  W 


Koe.Ea„  (195) 


Equations  (195)  and  (196)  result  in  vectors  of  length  2  that  are  obtained  from  the 
matrix  multiplication  of  the  2xN  matrix  of  equation  (190)  with  the  vectors  of  length  N 
from  equations  (109)  and  (110).  The  first  term  of  the  two  component  array  resulting 
from  equation  (195)  is 
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This  term  is  the  flux  exiting  the  right  hand  side  of  a  spatial  cell  in  a  single  direction.  It 
is  calculated  by  adding  the  contributions  from  fine  angular  quadrature  fluxes  entering 
the  left  hand  side  of  the  spatial  cell,  scattered  within  the  cell,  or  emitted  within  the  cell 
multiplied  by  their  respective  spatial  quadrature  coefficients.  Similarly  the  second  term 
of  equation  (195)  is  the  flux  exiting  the  left  edge 
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The  first  term  of  equation  (196)  is  cell  average  flux  streaming  rightward 
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It  is  calculated  by  adding  the  contributions  from  fine  quadrature  fluxes  entering  the  left 


hand  side  of  the  spatial  cell,  scattered  within  the  cell,  or  emitted  within  the  cell 
multiplied  by  their  respective  transport  coefficients.  Similarly  the  second  term 
resulting  from  equation  (196)  is  the  cell  average  flux  streaming  leftward 
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The  goal  is  to  transform  the  transport  equations,  (175)  and  (176),  which  have  N 
directions  and  are  not  practical  to  solve  into  equivalent  equations  with  only  two 
directions  that  are  practical  to  solve.  This  transformation  can  be  done  by  multiplying 
equations  (197)  through  (200)  by  an  appropriately  chosen  factor,  one  ,  obtained  from 
equations  (177)  through  (184).  As  an  example,  equation  (197)  and  (198)  can  be 
written  as 
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Note  that  these  two  equations  are  in  the  form  of  equations  (42)  and  (43)  for  two 
directions.  However,  the  transport  coefficients  are  not  calculated  from  the  coarse 
angular  quadrature.  They  are  calculated  from  the  fine  angular  quadrature  transport 
coefficients  and  fluxes.  The  two  direction  transport  effective  spatial  quadrature 
coefficients  are 
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Each  of  these  effective  transport  coefficients  can  be  thonght  of  as  a  weighted  snm  of  N- 


direction  spatial  qnadratnre  coefficients.  The  weights  are  the  terms  in  each  of  the 
brackets  of  the  above  eqnations  that  mnltiply  the  fine  angle  transport  coefficients.  Note 
that  the  snm  of  each  of  these  terms  over  its  respective  snmmation  index  is  one.  For 
instance  the  snm  of  the  weights  of  equation  (203)  are 
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For  this  reason  the  array  of  weights  for  the  positive  direction  fluxes  denoted  by 


(210) 


represents  a  distribution  of  flux  over  the  positive  direction  hemisphere.  The  array  of 
edge  flux  weights  for  the  positive  direction  for  cell  i  will  be  denoted  by  the  more 
convenient  notation 
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Additionally,  the  array  of  edge  flux  weights  for  the  negative  direction  for  cell  i  will  be 
denoted  by 
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The  flux  weight  array  for  the  cell  average  scattered  source  in  the  positive  direction  for 
cell  i  is 
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The  flux  weight  array  for  the  cell  average  scattered  source  in  the  negative  direction  for 
cell  i  is 
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The  array  of  flux  weights  for  the  cell  emission  source  in  the  positive  direction  for  cell  i 
is 
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The  array  of  flux  weights  for  the  cell  emission  source  in  the  negative  direction  for  cell  i 
is 
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If  equations  (199)  and  (200)  are  multiplied  by  the  appropriate  identity  formed  by 
rearranging  equations  (177)  through  (184),  as  was  done  for  the  edge  flux  the  two 
components  of  the  cell  average  flux  are 
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The  resulting  two-direction  effective  spatial  quadrature  coefficients  for  the  contributions 


of  edge  flux  scattered  source  and  emission  source  to  cell  average  source  are 
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The  weights  used  to  collapse  the  fine  angle  average  flux  spatial  quadrature  coefficients 
are  the  same  as  the  flux  weights  used  to  collapse  the  edge  flux  transport  coefficients. 
Equations  (203)  through  (224)  provide  a  mechanism  to  collapse  an  N  direction  angular 
quadrature  into  an  equivalent  two  direction  angular  quadrature  then  calculate  edge  flux 
and  average  flux  values  across  the  entire  spatial  domain.  If  the  edge  and  average  fluxes 
were,  known  these  ratios  could  be  correctly  calculated  and  the  number  of  particles 
exiting  a  cell  edge  can  be  calculated  directly.  Once  this  is  done,  the  particles  crossing  a 
cell  edge  in  the  direction  of  either  hemisphere  of  Figure  2  can  be  calculated.  It  remains 
to  apportion  this  flux  into  each  of  the  ordinates  of  the  fine  angle  quadrature. 

The  flux  and  source  weights  already  discussed  provide  the  mechanism  to  apportion  the 
coarse  quadrature  flux  found  from  solution  of  the  spatially  coupled  transport  system 
into  fine  directions.  A  is  used  to  denote  identify  fine-angle  fluxes  calculated  by 


apportioning  a  spatially-coupled  coarse-angle  flux.  For  instance,  represents  the 

apportioned  flux  in  ordinate  n  at  edge  This  flux  is  calculated  by  multiplying  the 

corresponding  coarse  angle  flux  with  its  appropriate  flux  weight  y/^  •_!  /  .  The 

’2  2 

array  of  positive  direction  flux  is  calculated  by  multiplying  the  positive  direction  coarse 
flux  by  the  array  of  positive  direction  flux  weights 
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An  angle  space  distribntion  iteration  is  shown  schematically  in  Fignre  3.  The  iteration 
begins  with  approximate  fine  angle  edge  flnx. 


Cell  Edge 


Cell  Edge 


Figure  3:  Apportioning  two  elements  into  N  elements  using  Iteration  Flux  Weights 

In  Fignre  3  there  are  fonr  initial  fine  flnxes  exiting  the  lower  left  hemispheres 
y/i  =  4,  ^^2  =  2,  ^3  =  1,  ^^4  =  3  .  These  flnxes  have  corresponding  weights  of 

fi  =  45/4  =  4'  These  weights  are  nsed  to  collapse  the  fine  angle  transport 
coefficients  to  two  directions  for  the  center  hemispheres.  These  two  direction 
coefficients  are  then  nsed  in  a  spatially  conpled  system  to  calculate  new  fluxes  that  have 
accounted  for  particle  scatter  across  the  problem  space.  In  this  example  this  results  in 
two  flux  values  exiting  the  top  hemispheres.  This  flux  is  then  apportioned  into  the 
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original  four  directions  using  flux  weights  as  shown  in  the  right  hemispheres.  This 
results  in  four  fluxes  of  =400,  =200,  ^3  =100,  ^4  =300  in  the  example  given. 

These  flux  solutions  may  not  be  correct,  but  they  are  of  the  right  order  of  magnitude 
because  they  account  for  infinite  particle  scatters  within  spatial  cells  and  across  the 
spatial  domain. 

The  use  of  approximate  fine  angle  fluxes  to  collapse  their  respective  spatial 
quadrature  coefficients  into  effective  spatial  quadrature  coefficients,  followed  by 
computation  of  transport  coefficients,  which  are  then  used  to  compute  spatially  coupled 
two  direction  fluxes,  which  are  apportioned  back  into  improved  estimates  of  fine  angle 
flux,  leads  to  an  iterative  scheme. 

The  vector  of  incoming  cell  edge  fluxes  for  I  spatial  cells  is  estimated.  Equation 
(112)  is  used  to  calculate  cell  outgoing  edge  fluxes  and  equation  (111)  is  used  to 
calculate  cell  average  fluxes.  These  calculated  fluxes  are  used  to  calculate  flux  weights 
with  equations  (211)  through  (216).  A  flux  solution  is  then  found  using  collapsed  the 
effective  transport  coefficients  for  the  flux  exiting  a  hemisphere  in  the  +//  direction  and 
the  flux  exiting  a  hemisphere  in  the  -/u  direction.  This  two  direction  flux  is  then 
apportioned  into  the  original  N  directions  for  a  flux  solution. 

The  flux  solution  arrived  at  with  estimated  flux,  calculating  approximate  flux 
weights  and  approximate  transport  coefficients  is  an  estimate  that  fully  accounts  for 
particle  scatters  across  the  spatial  domain  and  within  each  spatial  cell.  The  flux 
solution  is  not  dependent  on  numerous  source  iterations  to  account  for  the  scatters  to 
estimate  the  scattering  source.  It  fully  accounts  for  particle  scatters  without  iteration. 

This  straightforward  numerical  method  using  an  estimated  edge  flux  to  collapse 
cell  quadrature  coefficients,  use  these  to  convert  them  to  cell  transport  coefficients, 
calculate  a  two-direction  spatially-coupled  flux,  then  apportion  that  flux  back  into  the 
fine  ordinates  to  obtain  an  improved  estimate  of  edge  fluxes  was  the  basis  of  the 
iteration  scheme  developed  and  investigated.  This  transport  scheme  had  significant 
computational  efficiency  when  compared  with  SI.  Operationally  this  iteration  works  in 
the  following  way: 

-  calculate  fluxes  using  two-direction  effective  transport  coefficients, 

-  apportion  this  flux  into  fine  ordinates  as  an  initial  estimate  for  edge  fluxes. 
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-  calculate  updated  fine  angle  edge  fluxes, 

-  use  these  updated  edge  fluxes  to  calculate  weights 

-  use  the  weights  to  collapse  the  spatial  quadrature  coefficients  into  effective  two- 
ordinate  spatial  quadrature  coefficients 

-  use  these  to  form  the  effective  two-ordinate  transport  coefficients 

-  apportion  the  flux  calculated  in  this  way  into  fine  ordinates, 

-  calculate  a  better  estimate  of  edge  flux,  weights  and  transport  coefficients, 

-  proceed  until  convergence  tolerance  is  met. 

The  entire  scheme  amounts  to  iteration  on  flux  distribution  and  collapsed  transport 
coefficients.  The  algorithm  for  the  Angular-Spatial  Distribution  Iteration  (ASDI) 
method  using  a  step  characteristic  spatial  quadrature  is  shown  in  Algorithm  4. 


Initialize 

Obtain  an  initial  guess  for  cell  edge,  and  zeroth  moment  angular  flux 
Use  2  direction  quadrature  to  estimate  cell  coupling  coefficients 
Calculate  2  Direction  angular  flux 
Apportion  flux  into  to  N  directions 
Repeat 

update  cell  edge  flux 

solve  fine  angle  resolution  fluxes  within  spatial  cells 

use  flux  weights  to  generate  improved  coarse  angle  cell  coupling 

coefficients 

calculate  spatially  coupled  coarse  angle  edge  flux 

The  results  of  implementation  of  this  algorithm  in  Fortran  95  code  and 
comparison  with  SI  are  presented  and  discussed  in  the  next  chapter. 
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III.  sc  Experimental  Results 


I  set  out  to  design  a  transport  method  that  is:  robust,  efficient  (requires  few  iterations 
to  converge),  computationally  effective  (converges  rapidly  as  measured  by  compute 
time),  remains  so  when  confronted  with  material  discontinuities  and  is  readily 
parallelizable.  After  deriving  the  method,  developing  the  algorithms  and  implementing 
these  algorithms  in  FORTRAN  code  I  designed  experiments  to  test  4  of  the  5  desirable 
characteristics.  I  designed  experiments  to  test  the  methods  accuracy,  effectiveness  and 
computationally  efficiency  first  without  material  discontinuities  then  with  material 
discontinuities.  I  did  not  test  parallelizability  leaving  this  for  future  research. 

I  first  experimented  with  an  optically  thick  homogeneous  material.  I  varied  this 
homogeneous  material  from  absorber  (low  scattering  ratio)  to  scatterer  (scattering  ratio 
nearly  one)  in  order  to  test  ASDI  against  SI  for  two  types  of  problems:  one  in  which  SI 
converged  readily  and  one  for  which  it  converged  slowly. 

In  the  second  experiment  I  investigated  the  effect  of  periodic  material 
discontinuity  on  the  comparative  performance  of  ASDI.  To  do  this  I  chose  two 
materials  of  the  same  dimension.  I  placed  an  emission  source  in  one  material  and  no 
source  in  the  other  material.  Both  materials  were  I  MFP  wide.  The  emitter  had  a 
scattering  ratio  of  one.  The  non-emitter’s  scattering  ratio  varied  from  0.0  to  I.O  as  in 
the  homogeneous  material.  This  two  material  pattern  was  repeated  10  times  creating  a 
periodic  material  discontinuity. 

In  addition  to  scattering  ratio  I  further  investigated  the  effect  of  spatial 
refinement,  angular  refinement  and  convergence  criterion  on  robustness,  effectiveness, 
and  efficiency  for  both  the  homogeneous  and  periodic  problems.  I  did  this  by  setting 
the  scattering  ratio  to  0.98  in  the  homogeneous  material  of  problem  I  and  the  non¬ 
emitter  of  problem  2  instead  of  allowing  this  parameter  to  vary.  The  parameters  I  then 
varied  in  turn  were  spatial  refinement,  angular  refinement  and  convergence  tolerance. 

I  tested  ASDFs  accuracy  by  checking  its  solution  against  a  benchmark  to  ensure 
it  met  the  convergence  tolerance  required.  The  first  benchmark  I  used  was  an 
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unaccelerated  SI  allowed  to  converge  to  the  same  convergence  tolerance  as  ASDI.  This 
provided  a  straightforward  and  reliable  benchmark  but  only  for  test  problems  that  were 
absorptive  in  which  SI  did  not  falsely  converge.  For  those  problems  in  which  SI  did 
falsely  converge  I  took  advantage  of  Si’s  ability  to  recognize  a  correct  answer  instead  of 
its  ability  to  calculate  a  correct  answer.  First  an  initial  guess  of  zero  was  used  and  a 
conventional  SI  solution  was  found.  This  solution  could  be  off  by  as  much  as  four 
orders  of  magnitude  because  of  false  convergence.  In  this  case  the  flux  solution  was 
discarded  but  the  iterations  required  to  obtain  the  solution  were  recorded.  I  then  input 
ASDI  solution  to  SI  as  an  initial  guess  and  allowed  SI  to  iterate  the  same  number  of 
times  it  used  to  obtain  its  falsely  converged  solution.  My  assumption  was  that  SI  would 
recognize  a  fixed  point.  I  allowed  SI  to  iterate  the  same  number  of  times  it  took  to 
obtain  its  falsely  converged  solution  to  ensure  that  it  had  enough  iteration  to  drift  away 
from  the  ASDI  solution  if  the  initial  value  (i.e.  the  ASDI  solution)  was  not  a  fixed 
point.  For  instance  if  SI  required  100000  iterations  to  converge  on  an  answer  I  recorded 
this  but  discarded  its  flux  solution.  I  then  passed  SI  the  ASDI  solution  and  allowed  it 
to  iterate  100000  times  without  checking  it  for  convergence.  After  completing  these 
100000  iterations  the  new  SI  (SI_ASDI)  solution  was  compared  with  the  ASDI  solution 
using  (  ASDI,SI_ASDI)  .  This  symmetric  relative  difference  was  then  compared  with 

the  desired  convergence  tolerance.  If  it  was  tighter  than  this  desired  convergence 
tolerance  it  demonstrated  reliable  accuracy.  After  ensuring  accuracy  effectiveness  was 
tested,  I  examined  computational  efficiency.  I  tested  this  by  comparing  ASDI  and  SI 
compute  times. 


Problem  1  Optically  Thick  Homogeneous  Material 


The  first  test  problem  studied  is  a  100  cm  thick  homogeneous  slab  with  a  vacuum 
boundary  on  the  right,  a  symmetry  boundary  on  the  left,  a  uniform  isotropic  source. 
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source,  and  a  total  cross  section  of  1.0  cm'^.  A  diagram  of  this  problem  is  shown  in 
Figure  4. 


=  0.0 


Figure  4  Problem  1 

In  the  study  of  this  problem  and  other  problems  the  terms  refinement  factor  and 
scattering  ratio  fraction  are  introduced.  Refinement  factor  is  an  integer  that  describes 
the  number  of  partitions  chosen  for  the  spatial  domain  in  a  mesh.  For  instance  the 
spatial  domain  of  problem  1  was  100  cm.  Refinement  factors  used  were  {1,  2,  4,  8,  10, 
16,  32,  64,  100,128}  designating  the  number  of  cell  partitions.  These  refinement  factors 

correspond  with  cell  widths  of  {  100,50,25,12.5,10,^,  Problem  1. 

Scattering  ratio  fraction  is  a  real  number  used  as  a  coefficient  that  multiplies  a  base 
scattering  ratio  defined  by  the  user.  Scattering  ratio  fraction  is  used  order  to  vary 
problem  scattering  ratio  over  a  range  of  interest  for  the  numerical  experiment.  For 
instance,  in  problem  1  the  base  scattering  ratio  was  1.0.  The  scattering  ratio  fraction 
varied  from  0.0  through  1.0.  These  factors  were  used  to  vary  the  parameters  under 
study  across  a  full  range  of  interest. 

The  ASDI  method  provides  a  useful  transport  tool  only  if  it  provides  reliable  and 
usably  accurate  numerical  solutions.  The  first  experiment  tested  the  ASDI  solution 
robustness  as  scattering  ratio  varied.  Figure  5  displays  two  plots.  These  are  the 
symmetric  relative  difference  versus  scattering  ratio  of  the  ASDI  solution  and  the  SI 
solution  (ASDI_SI)  and  of  the  ASDI  solution  and  the  SI  solution  after  being  given  the 
ASDI  solution  (ASDI_SI-ASDI).  The  convergence  tolerance  chosen  was  10^®  for  both 


3-3 


SI  and  ASDI.  At  a  scattering  ratio  of  0.6  the  ASDI  and  SI  solntions  differ  by  more 
than  the  chosen  convergence  tolerance.  This  difference  increases  as  the  scattering  ratio 
increases.  Both  methods  agree  within  the  convergence  tolerance  for  scattering  ratios 
below  0.6.  These  plots  show  that  SI  produces  unreliable  answers  at  scattering  ratios  as 
low  as  0.6  but  ASDI  continues  to  provide  the  accurate  answers  are  validated  by  STs 
fixed  point  recognition  for  all  scattering  ratios. 


Figure  5:  Symmetric  Relative  Difference  between  the  ASDI  solution  and  the  SI 
solution  as  scattering  ratio  varies.  Angular  quadrature  is  DE-8,  refinement  is  50 
(crAx  =  2MFPs),  convergence  tolerance  is  10 


Figure  6  shows  iteration  count  versus  scattering  ratio  (c)  with  an  angular 
refinement  of  8,  a  mesh  width  of  2  MFPs  and  a  convergence  tolerance  of  10~®.  SI 
iteration  count  increases  with  c,  climbing  steeply  as  c  approaches  unity.  This 
demonstrates  the  classic  weakness  of  SI  for  diffusive  problems.  The  ASDI  solution 
converges  in  two  or  three  iterations  regardless  of  scattering  ratio  demonstrating  its 
usefulness  especially  for  diffusive  problems. 
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Figure  6  Plot  of  iteration  count  versus  scattering  ratio.  Angular  quadrature  is  DE-8, 


refinement  is  50  (Ax  =  2MFPs  ),  convergence  tolerance  is  10 


The  low  iteration  count  of  the  ASDI  method  displayed  in  Figure  6  is 
encouraging.  It  is  not  a  practical  improvement  over  SI  unless  computational  cost  is  also 
reduced.  Figure  7  displays  compute  time  versus  scattering  ratio.  The  amount  of 
compute  time  needed  by  the  ASDI  method  is  insensitive  to  the  scattering  ratio  and 
compares  favorably  to  SI  even  when  SI  converges  rapidly.  This  demonstrates  that  we 
have  not  paid  excessive  computational  efficiency  costs  to  obtain  the  method’s  desirable 
robustness  and  effectiveness. 
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Figure  7  Plot  of  compute  time  (seconds)  versus  scattering  ratio.  Angular  quadrature 
is  DE-8,  refinement  is  50  {a Ax  =  2MFPs),  convergence  tolerance  is  10~®. 


The  results  shown  in  Fignre  5  throngh  Fignre  7  demonstrate  the  effect  of  scattering 
ratio  on  the  ASDI  and  SI  methods.  These  results  were  obtained  with  a  fixed  refinement 
of  50  mesh  cells  {Ax  =  2.0cm),  convergence  tolerance  fixed  atl0~®,  and  angular 
quadrature  fixed  at  n  =  8  .  Spatial  refinement,  angular  refinement,  and  convergence 
tolerance  are  also  parameters  that  are  expected  to  have  an  impact  on  robustness, 
effectiveness,  and  efficiency  of  a  method.  The  next  series  of  plots  investigates  the  effect 
of  changing  these  parameters.  These  three  parameters  were  studied  with  a  scattering 
ratio  of  0.98  making  the  homogeneous  material  a  good  scatterer.  This  scattering  ratio 
was  chosen  to  test  the  performance  of  ADSI  in  a  diffusive  problem. 

Figure  8  shows  the  relative  difference  between  SI  and  ASDI  solutions  as  the 
spatial  mesh  is  refined.  The  difference  between  the  answers  is  greater  than  the 
convergence  tolerance  because  with  a  scattering  ratio  of  0.98  SI  suffers  from  false 
convergence.  If  SI  is  fed  the  ASDI  solution  it  again  recognizes  this  solution  as  a  fixed 
point  for  all  spatial  refinements. 
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Figure  8  Plot  of  Symmetric  Relative  Difference  £srd  between  ASDI  solution  and  SI 
solutions  for  as  cell  mesh  is  refined.  Angular  quadrature  is  DE-8,  scattering  ratio  is  1.0, 

convergence  tolerance  is  10  ® . 


SI  recognition  of  the  ASDI  solntion  as  a  fixed  point  for  all  spatial  refinements 
provides  convincing  evidence  of  the  accnracy  of  its  solntion.  I  was  able  to  farther  test 
this  accnracy  because  by  computing  and  comparing  the  ASDI  convergence  rate  with  the 
known  convergence  rate  of  SC.  An  analytic  solution  was  available  for  problem  1.  This 
was  done  for  a  DE-4  angular  quadrature  and  a  scattering  ratio  of  1.0  with  spatial 
refinement  varying  from  1000  (cell  width  =  1.25xl0~^cm)  through  10  (cell 
width=10cm).  ASDI  solutions  are  not  expected  to  be  the  same  as  the  analytic  solution 
because  ASDI  is  a  numerical  approximation  of  the  analytic  solution  dependent  on  mesh 
refinement.  However,  because  SC  is  known  to  be  second  order  convergent  in  space  the 
method  can  be  checked  against  the  analytic  solution  to  determine  if  the  order  of 
convergence  is  2.  Figure  9  shows  symmetric  relative  difference  {ssRD )  between  the 
ASDI  solution  and  an  exact  solution  plotted  against  the  spatial  mesh  refinement  factor 
[Rp)  on  a  Log-Log  graph  for  a  scattering  ratio  of  1.0.  Figure  10  shows  this  for  a 
scattering  ratio  of  0.9.  egRD  appears  as  a  straight  line  on  both  plots  whose  slope  is  two. 
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A  line  with  a  slope  of  two,  which  is  the  convergence  order  of  SC,  has  been  overlaid  on 
the  data  plot.  This  agreement  of  ASDI’s  convergence  rate  with  the  known  SC 
convergence  rate  for  a  scattering  ratio  of  1.0  and  0.9  combined  with  SI  recognition  of 
ASDI  solntions  as  fixed  points  demonstrate  the  method’s  robnstness. 


Figure  9:  Plot  of  Symmetric  Relative  Difference  between  the  ASDI  solution  and 
an  analytic  solution  as  cell  mesh  is  refined  between  10  MFPs  and  0.1  MFPs.  Angular 
quadratme  is  DE-8,  scattering  ratio  is  1.0,  convergence  tolerance  is  10~® . 
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Figure  10:  Plot  of  Symmetric  Relative  Difference  egjuj  between  the  ASDI  solution  and 
an  analytic  solution  1  as  cell  mesh  is  refined  between  10  MFPs  and  0.1  MFPs.  Angular 
quadrature  is  DE-8,  scattering  ratio  is  0.9,  convergence  tolerance  is  10~® . 

Figure  11  shows  the  iteration  count  of  the  ASDI  and  SI  methods  as  the  mesh  is 
refined.  SI  iteration  remains  flat,  at  500  iterations,  as  the  mesh  is  refined  and  ASDI 
iteration  count  increases  from  4  iterations  for  a  coarse  mesh  to  17  iterations  for  a  fine 
mesh  indicating  ASDI  is  an  effective  method.  SI  appears  insensitive  to  the  cell 
refinement.  I  drew  no  conclusions  about  this  SI  characteristic  from  this  because  of  its 
false  convergence. 
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Figure  11  Plot  of  iteration  count  as  cell  mesh  is  refined.  Angular  quadrature  is  DE-8, 
scattering  ratio  is  1.0,  convergence  tolerance  is  10~® . 

Figure  12  shows  the  compute  times  required  for  the  SI  and  ASDI  solutions  as  the 
spatial  mesh  is  refined.  The  spatial  refinement  increases  compute  time  for  both 
methods.  ASDI  requires  less  compute  time  than  SI  for  the  scattering  ratio  chosen.  This 
indicates  that  ASDI  will  not  require  disproportionately  more  time  than  SI  regardless  of 
mesh  size.  The  true  value  of  ASDI  is  that  it  converges  to  the  correct  solution  without 
requiring  large  compute  time,  something  that  SI  can  not  do. 
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Figure  12  Plot  of  compute  time  required  for  the  ASDI  and  SI  solutions  as  cell  mesh  is 
refined.  Angular  quadrature  is  DE-8,  scattering  ratio  is  1.0,  convergence  tolerance  is 

10"®. 

The  convergence  tolerance  nsed  to  test  previons  problems  was  10"®  on  angnlar 
flnx.  This  is  a  fairly  tight  convergence  tolerance  and  it  is  snitable  for  most  engineering 
applications.  However  some  applications  may  reqnire  tighter  tolerances.  The  effect  of 
tightening  convergence  tolerance  from  10"® to  10"^^  on  symmetric  relative  difference, 
iteration  connt,  and  compnte  time  is  shown  next.  Fignre  13  shows  that  the  accnracy  of 
the  converged  SI  solntion  does  not  meet  the  accnracy  reqnired  by  the  specified 
convergence  tolerance  regardless  of  how  tight  that  tolerance  is,  whereas  the  ASDI 
method  continnes  to  provide  reliably  accnrate  solntions  for  any  tolerance  withont  ronnd 
off  error  or  instability.  Examination  of  Fignre  14  and  Fignre  15  reveals  that  this 
accnracy  is  achieved  with  modest  increase  in  iterations  reqnired  or  compnte  time. 
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Figure  13  Plot  of  symmetric  relative  difference  versus  convergence  tolerance.  Angular 
quadrature  is  DE-8,  scattering  ratio  is  1.0,  cell  size  is  2  MFPs. 
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Figure  14  Plot  of  iterations  versus  convergence  tolerance  for  ASDI  and  SI.  Angular 
quadrature  is  DE-8,  scattering  ratio  is  1.0,  cell  size  is  2  MFPs. 
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Figure  15  Plot  of  compute  time  versus  convergence  tolerance  for  ASDI  and  SI. 

Angular  quadrature  is  DE-8,  scattering  ratio  is  1.0,  cell  size  is  2  MFPs. 

The  effect  of  increasing  angnlar  refinement  is  shown  in  Fignre  16  throngh  Fignre 
18  shows  that  SI  remains  inaccnrate  and  has  increased  compntational  cost  as  the 
angnlar  qnadratnre  is  refined.  ASDI  continnes  to  provide  reliably  accnrate  solntions 
withont  an  increase  in  iteration  connt  or  compnte  time  beyond  an  angnlar  qnadratnre  of 
2.  ASDI  achieves  a  quick  solution  for  an  angular  quadrature  of  2  because  its  global 
spatial  solution  was  designed  with  this  angular  quadrature.  Such  a  coarse  angular 
quadrature  is  not  likely  to  meet  most  engineering  needs.  However  ASDFs  rapid 
convergence  for  even  fine  angular  quadratures  indicates  usefulness  for  engineering 
problems  requiring  higher  angular  refinement. 
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Figure  16  Plot  of  symmetric  relative  difference  versus  angular  quadrature. 
Convergence  tolerance  is  10  ®  scattering  ratio  is  1.0,  cell  size  is  2  MFPs. 


Figure  17  Plot  of  iterations  versus  angular  quadrature  for  ASDI  and  SI.  Convergence 
tolerance  is  10~®  scattering  ratio  is  1.0,  cell  size  is  2  MFPs. 
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Figure  18  Plot  of  compute  time  versus  angular  quadrature  for  ASDI  and  SI. 

Convergence  tolerance  is  10  ®  scattering  ratio  is  1.0,  cell  size  is  2  MFPs. 

This  section  demonstrated  the  strength  of  the  ASDI  method  for  a  homogeneons 
material.  The  most  important  characteristic  examined  was  accnracy  for  each  of  the 
parameters  stndied.  Each  experiment  showed  that  ASDI  was  reliably  accnrate  within 
the  convergence  tolerance  reqnired  whereas  SI  was  not.  This  reliable  accnracy  across  all 
parameters  stndied  demonstrates  that  the  method  is  robnst  for  this  homogeneons 
problem.  The  next  section  examines  the  method’s  performance  when  applied  to  a 
heterogeneons  material. 


Periodic:  Two  Regions  Repeated  10  Times 


The  second  test  problem  investigated  was  slab  with  a  symmetry  bonndary  on  the  left 
side  and  a  vacnnm  bonndary  on  the  right  side.  Two  materials  of  1  cm  cell  width  are 
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placed  side  by  side.  This  two  material  pattern  is  repeated  10  times  for  a  total  length  of 
20  cm.  Material  I  had  a  total  cross  section  of  1.0  cm'^ ,  a  scattering  ratio  of  1.0,  and  a 

-3 

uniform  source  of  1.0  cm"  .  Its  material  properties  remain  fixed  for  all  experiments. 
Material  II  has  a  total  cross  section  of  1.0  cm'^ ,  a  baseline  scattering  ratio  of  1.0,  and  no 
source.  Material  IPs  parameters  were  varied  during  the  experiments.  A  diagram  of 
these  two  materials  is  shown  in  Figure  19. 
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Figure  19  Problem  2 

In  the  next  set  of  experiments  I  tested  the  method  for  accnracy,  effectiveness,  and 
efficiency  for  varying  scattering  ratio  fraction,  refinement  factor,  convergence  tolerance 
and  angnlar  refinement  as  was  done  in  problem  1. 

The  first  set  of  experiments  for  problem  2  tested  the  ASDI  method’s  accnracy, 
effectiveness  and  efficiency  versns  scattering  ratio  fraction  in  Material  11.  Scattering 
ratio  fraction  was  varied  from  0.0  to  1.0.  Angnlar  refinement  for  this  experiment  was  8, 
refinement  factor  was  1  and  convergence  tolerance  was  10~® . 

Accnracy  is  the  first  parameter  presented.  Fignre  20  displays  the  symmetric 
relative  difference  between  ASDI  and  SI  (ASDI,SI_ASDl))  and  between  ASDI  and 
SI  given  the  ASDI  solntion  (ASDI,SI_ASDl)).  The  plots  show  that  ASDI  provides 
reliably  accnrate  solntions  for  all  scattering  ratios  while  SI  provides  reliably  accnrate 
solntions  only  throngh  scattering  ratios  of  0.6. 
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Figure  20:  Symmetric  Relative  Difference  between  the  ASDI  solution  and  the  SI 
solution  as  scattering  ratios  vary.  Angular  quadrature  is  DE-8,  convergence  tolerance 

is  10  ® ,  refinement  is  50. 

The  effectiveness  of  the  method  is  tested  next.  Fignre  21  shows  the  method 
iteration  connt  versns  scattering  ratio  fraction  for  material  11.  This  plot  shows  that 
ASDI  converges  in  fonr  or  five  iterations  regardless  of  scattering  ratio  demonstrating  its 
effectiveness  for  diffnsive  heterogeneons  materials.  Iteration  connt  for  SI  increases  with 
scattering  ratio.  This  connt  climbs  steeply  as  scattering  ratio  approaches  1.0. 
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Figure  21:  Plot  of  iteration  count  versus  scattering  ratio  fraction.  Angular  quadrature 
is  DE-8,  convergence  tolerance  is  10  ® ,  refinement  is  50. 
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Efficiency  of  the  method  is  tested  next.  Fignre  22  shows  compnte  time  versns 
scattering  ratio  fraction  for  material  11.  This  fignre  shows  that  ASDI  compnte  time  is 
less  than  0.1  seconds  for  all  scattering  ratios  tested.  Comparison  of  the  ASDI  and  SI 
compnte  times  demonstrate  that  ASDI  is  more  compntationally  efficient  than  SI  even 
when  material  II  is  a  strong  absorber  and  SI  converges  rapidly.  The  fignre  also  shows 
that  in  materials  that  are  strong  scatterers  SI  compnte  time  climbs  rapidly. 
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Figure  22:  Plot  of  compute  time  (seconds)  versus  scattering  ratio.  Angular  quadrature 


is  DE-8,  convergence  tolerance  is  10  ,  refinement  is  50. 


The  next  set  of  experiments  demonstrate  method  accnracy,  effectiveness  and 
efficiency  while  refinement  factor  varies  for  material  11.  Scattering  ratio  for  material  II 
was  fixed  at  1.0.  This  tested  ASDI  in  a  highly  diffnsive  material.  A  material 
discontinnity  occnrs  even  thongh  cross  sections  and  scattering  ratios  are  the  same 
becanse  intrinsic  sonrces  are  emitted  only  in  material  I.  The  refinement  factor  for 
material  II  was  varied  between  1  and  66. 

Accnracy  was  tested  first.  Fignre  23  displays  the  symmetric  relative  difference 
between  ASDI  and  SI  (ASDI,SI_ASDl))  and  between  ASDI  and  SI  given  the  ASDI 
solntion  (£’gj^j)(ASDI,SI_ASDl)).  The  plots  show  that  ASDI  provides  reliably  accnrate 
solntions  for  all  refinement  factors  bnt  SI  does  not. 
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Figure  23:  Plot  of  Symmetric  Relative  Difference  between  ASDI  solution  and  SI 
solutions  as  cell  mesh  is  refined.  Angular  quadrature  is  DE-8,  scattering  ratios  are 
baseline  values,  convergence  tolerance  is  10~® . 


Effectiveness  was  next  tested.  Fignre  24  shows  that  ASDI  converges  on  a  solntion  in 
less  than  8  iterations  for  a  refinement  factor  of  1  growing  slightly  as  refinement  factor  is 
increased  then  remaining  constant  for  refinement  factors  above  16.  SI  reqnires  many 
more  iterations  than  ASDI  and  provides  nnreliable  solntions  for  these  refinement  factors 
since  this  is  a  diffnsive  problem. 
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Figure  24  Plot  of  iteration  count  as  cell  mesh  is  refined.  Angular  quadratme  is  DE-8, 
scattering  ratios  are  baseline  values,  convergence  tolerance  is  10  ® . 

Computational  efficiency  was  next  tested.  Figure  25  demonstrates  that  ASDI 
requires  less  compute  time  than  SI  regardless  of  mesh  size  for  the  baseline  scattering 
ratios  selected.  However,  as  spatial  mesh  is  refined  the  compute  time  of  ASDI  is  nearly 
the  same  as  the  compute  time  of  SI.  The  relatively  comparable  compute  times  of  SI 
and  ASDI  do  not  indicate  that  SI  is  just  as  efficient  as  ASDI  for  these  fine  spatial 
meshes  because  SI  does  not  provide  reliably  accurate  answers.  ASDI  does  not  require 
excessive  compute  time  at  all  refinement  factors. 
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Figure  25:  Plot  of  compute  time  as  cell  mesh  is  refined.  Angular  quadrature  is  DE-8, 
scattering  ratios  are  baseline  values,  convergence  tolerance  is  10  ® . 


Accuracy,  effectiveness  and  compntational  efficiency  were  fnrther  tested  with 
fixed  scattering  ratio  and  refinement  factor  for  varying  convergence  tolerance  and 
angnlar  refinement.  These  resnlts  are  not  shown  here.  As  in  problem  1  these  resnlts 
indicate  that  ASDI  remains  reliably  accnrate,  effective  and  efficient  as  these  parameters 
vary. 

These  experiments  demonstrate  the  strength  of  the  ASDI  method  for  a  2  region 
periodic  material.  The  experiments  demonstrate  that  ASDI  was  always  reliably 
accnrate  and  that  the  parameters  with  the  greatest  impact  on  iteration  connt  and 
compnte  time  were  scattering  ratio  and  refinement  factor.  The  experiments  presented 
kept  cross  sections  constant  bnt  varied  spatial  refinement.  In  order  to  fnlly  stress 
material  discontinnity  in  this  periodic  material  I  next  examined  the  impact  of  varying 
both  scattering  ratio  and  cross  section  simnltaneonsly  in  material  II.  I  varied  scattering 
ratio  fraction  from  0.0  to  1.0  and  cross  section  from  10~^cm~^  to  lO^^cm”^.  I  examined 
the  impact  of  varying  these  two  parameters  on  accnracy,  effectiveness  and  efficiency  at 
a  refinement  factor  of  1. 

A  3D  plot  of  ASDI  the  symmetric  relative  difference  between  ASDI  and  SI 
( ^SRD  (ASDI,SI_ASDI)  )  and  between  ASDI  and  SI  given  the  ASDI  solntion 
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( ^SRD  ("^SDI,SI_ASDl)  )is  shown  in  Fignre  26.  This  plot  shows  that  ASDI  is  reliably 
accnrate  across  the  wide  range  of  material  properties  tested 
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Periodic  SC  Refinement  =1 


10-7 

Figure  26:  3D  Plot  of  ASDI  symmetric  relative  difference  between  ASDI  and  SI 
( ^SRD  ( ASDl,Sl_ASDl) )  and  between  ASDI  and  SI  given  the  ASDI  solution 
( ^SRD  ( ASDl,Sl_ASDl) )  as  scattering  ratio  and  cross  section  vary  .  Angular  quadrature 
is  DE-8,  convergence  tolerance  is  10  ® ,  refinement  factor  is  1. 


A  3D  plot  of  ASDI  iteration  count  versus  cross  section  and  scattering  ratio  is 
shown  in  Figure  27. 
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Periodic  SC  Refinement  =1 


Figure  27:  3D  Plot  of  ASDI  iteration  count  as  scattering  ratio  and  cross  section  vary 


Angular  quadrature  is  DE-8,  convergence  tolerance  is  10  ® ,  refinement  factor  is  1. 


This  figure  demonstrates  that  ASDI  is  not  sensitive  to  either  scattering  ratio  or  cross 
section  taking  at  most  7  iterations  to  converge  regardless  of  the  sharpness  of  material 
discontinnity.  Again,  a  low  iteration  connt  is  not  usefnl  if  compntational  cost  is 
prohibitive.  A  comparison  of  the  time  reqnired  for  ASDI  and  SI  to  solve  the  problem  is 
shown  in  Fignre  28.  The  ratio  of  the  ASDI  compnte  time  to  SI  compnte  time  is  shown 
in  the  fignre,  if  this  ratio  is  less  than  1.0  then  ASDI  takes  less  time  than  SI.  This  is  the 
case  for  almost  all  combinations  of  scattering  ratio  and  cross  section.  SI  takes  less 
compnte  time  than  ASDI  when  scattering  ratios  that  are  not  diffnsive.  Even  in  these 
problems  where  SI  is  computationally  efficient  ASDI  compute  times  is  of  the  same  order 
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of  magnitude  but  more  importautly  provides  reliably  accurate  auswers  across  a  broader 
rauge  of  parameters.  SI  provides  uureliable  auswers  because  it  couverges  falsely  eveu 
wheu  material  II  is  a  stroug  absorber.  The  false  couvergeuce  of  SI  occurs  because 
material  I  is  diffusive. 

Periodic  SC  Refinement  =1 


Figure  28:  3D  Plot  of  the  ratio  ASDI  compute  time  to  SI  Compute  time  as  scattering 
ratio  and  cross  section  vary  .  Angular  quadrature  is  DE-8,  convergence  tolerance  is 

10  ® ,  refinement  factor  is  1. 
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These  experiments  demonstrate  that  ASDI  is  reliably  accurate  for  all  the 
parameters  studied  and  therefore  robust.  The  method  dramatically  reduces  iteration 
count  and  required  computing  time  for  diffusive  problems  compared  with  SI. 
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IV  Solution  of  the  Discrete  Ordinates  {S„)  Transport  Equations  with 

EC 


This  chapter  introduces  cell  transport  coefficients  for  the  EC  spatial 
quadrature.  The  derivation  of  this  spatial  quadrature  was  done  by  Mathews, 

Minor,  and  Sjoden  in  1993.  I  adopt  their  notation  for  limits  of  integration  and 
leftward  or  rightward  directed  fluxes  and  flux  moments.  With  this  notation  the 
cell  average  source  is 

(227) 


(228) 


where  Av^  is  the  width  of  a  spatial  cell  whose  left  edge  is  x._i  and  whose  right 

'  2 


edge  is  the  x  j  .  The  average  scattered  source  in  cell  i  is 
/h — 

2 


n4 


The  average  intrinsic  particle  emission  in  cell  i  is 


(229) 


(230) 


The  SC  spatial  quadrature  assumed  that  the  both  source  distributions  were 
constant.  However,  the  EC  spatial  quadrature  assumes  that  the  distribution  of 
scattered  source  is  exponential.  The  distribution  of  emission  source  is  assumed  to 
be  constant  for  this  discussion.  An  alternative  to  a  constant  emission  source  is  to 
assume  it  has  the  same  exponential  distribution  as  the  scattered  particle  source 
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distribution.  If  this  is  done  the  averaging  of  the  emission  source  is  handled  in  the 
same  way  as  the  scattered  source.  Extension  for  an  exponentially  distributed 
emission  source  is  obvious. 


EC  Transport  coefficients 


The  scattered  source  distribution  for  EC  is  exponential 


(231) 

The  parameters  j  and  j  are  chosen  to  match  available  information  about 
the  source,  typically  its  average,  equation  (38)  and  first  spatial  moment 
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where  Pi(x)  is  a  shifted  Legendre  polynomial: 


(232) 


Pi(x)  =  2 


1, 


(233) 


defined  in  the  interval  to  In  practice  A  ^  ,5'„  are  not  calculated 

from  equations  (38)  and  (232)  but  are  accumulated  from  angular  flux  moments 
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EC  transport  equations  for  cell  edge  flux  are 


^nj^i  Pn  i  ^/f  (I  I  \ 

r„,,i  =r„,,._ie '  ■  X(Km|-A.,). 
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Mn  <  0  (237) 


See  page  31  of  the  Mathews,  Sjoden  and  Minor  work  (reference  7)  where  these 
equations  are  clearly  derived.  Optical  path  length,  ^ ,  was  given  previously  as 
equation  (56).  The  transport  equations  for  cell  average  flux  are 
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The  transport  equations  for  cell  first  spatial  flux  moments  are 
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(240) 


+ 


\Mn\  L 


(241) 


Examination  of  these  equations  reveals  they  are  in  the  form  of  equations  (42) 
through  (45).  The  transport  coefficients  are  obvious.  The  first  transport 
coefficient  from  the  first  term  of  equation  (236)  or  (237)  relates  outgoing  flux  to 
incoming  flux.  It  is 
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(242) 
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The  second  transport  coefficient  from  the  second  term  of  eqnation  (236)  or  (237) 
describes  the  contribntion  of  within  cell  scattering  to  the  cell  ontgoing  flnx.  It  is 


The  third  transport  coefficient  from  the  first  term  of  eqnations  (238)  or  (239) 
describes  contribntion  of  flnx  entering  a  cell  edge  to  cell  average  flnx.  It  is 


(243) 


(244) 

The  fonrth  transport  coefficient  from  the  second  term  of  eqnations  (238)  or  (239) 
describes  contribntion  of  cell  scattering  to  cell  average  flnx.  It  is 


(246) 

The  fifth  and  sixth  transport  coefficients  depend  on  whether  a  particle  enters  a 
cell  right  or  left  edge.  These  coefficients  captnre  flnx  spatial  first  moment 
information  the  negative  sign  on  equation  (241)  reflect  that  a  flux  gradient 
appears  different  to  particles  streaming  in  ordinates  that  have  equal  direction 
cosines  but  in  opposing  directions.  The  fifth  transport  coefficient,  from  the  first 
term  of  equations  (240)  or  (241),  describes  contribution  of  flux  entering  an  edge 
to  cell  first  moment  flux.  It  is 
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The  sixth  transport  coefficient,  from  the  second  term  of  equations  (240)  or  (241), 
describes  contribution  of  cell  scattering  to  cell  first  moment  flux.  It  is 
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(247) 


Comparison  of  equations  (242)  through  equation  (245)  shows  that  the  difference 
between  SC  transport  coefficients  and  EC  transport  coefficients  is  that  Ani  =  o  for 
SC  because  its  spatial  first  moment  is  defined  to  be  zero. 

The  cell  EC  transport  equations  may  be  written  in  vector  notation  as 
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where 


y^out^  y^in^  '^^OSj  (Pi)  +^0£';  ^A  ’  (2^®) 

V^A^  ^^Aii  V^i„i  +K^5.  ([3;)  Sai  +K^£.  Eai  ,  (249) 

=Kx/,  Win,  +K^.  54.  ,  (250) 

SA,=^^WAi^  (251) 

(252) 


A. 
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PN,i 


where  each  of  the  coefficient  matrices  is  dependent  on  the  associated  element  of 
[3;.  For  the  same  reason  that  assembly  of  a  fnlly  spatially  and  directionally 
conpled  system  of  equations  was  impractical  to  solve  for  SC  it  is  even  more 
impractical  for  EC.  This  is  because  of  the  additional  first  moment  equations 
(250)  and  the  inherent  nonlinearity  in  (7:29,30  ). 

However,  if  the  incoming  edge  flux  is  known  can  be  found  by  iteration.  The 
logic  of  this  iteration  is  shown  in  Algorithm  5. 
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Initialize 

estimate  edge  flux 

calculate  cell  coupling  coefficients  for 


Do 


calculate 


A  ’  r 
V  ^  X  ^ 


calculate  K 


V  y 


calculate  y/^f 
calculate  |  Vm  j 


calculate  S 


calculate  S 


calculate  new  ,5*^ 


End  do  when  SsRoi^f  converged 


Equations  (89)  through  (252)  are  used  in  the  beta  iteration  algorithm. 


Once 


Pi  is  known,  edge  flux  is  also  found  by  iteration.  Outgoing  edge  flux  is 
calculated  from  equation  (89)  and  y/ is  calculated  from  (90).  This  is  done  for 
each  cell  and  then  boundary  conditions  are  applied.  The  outgoing  edge  flux 
calculated  for  a  cell  is  used  as  the  incoming  edge  flux  for  an  adjacent  cell.  The 
logic  of  this  iteration  is  shown  in  Algorithm  6. 


4-6 


initialize  flux 


P  ^converged 


-(0)  -(Z,) 

y^n  =Vin 


Do  all  cells 


calculate  for  all  cells 


apply  boundary  conditions 
update  from  adjacent  cell  outflow 


End  do  when  SsRoiy'^^  current  P 


Analysis  of  Algorithm  5  and  Algorithm  6  reveals  that  these  two  algorithms 
are  linked  through  P  •  and  edge  flux  .  Changing  P  leads  to  new  transport 
coefficients  which  lead  to  new  flux  w  ■  .  Updated  y/-^  leads  to  new  estimates  for 

I  Hi  Itlj 

P  .  The  interdependence  of  these  two  algorithms  suggests  the  following  iteration 


scheme  to  couple  outgoing  edge  flux  iteration  with  beta  iteration. 
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Initialize 


-(0) 

estimate  initial  flux  \|/^-^ 


-(0)  -(*) 

initialize  flux  If/ 


Do  all  cells 


do) 

estimate  p 


calculate  cell  coupling  coefficients  for 


k(p'"> 


calculate 


^ ’K  (p)"]  for  all  cells 


apply  boundary  conditions 


— (/” +l) 

update  from  adjacent 

cell  outflow 


calculate  P^  ^  f  *^[4  ^ 


calculate  K  f  P^  ^ 


f~(f+l)  -(/■) 

Enddowhen  SsRD  Win  ^W in 

V 

for  current  P 


calculate 


-(/■+1)  /-(Z,) 
calculate  y/^^  M 


calculate  S 


calculate  S 


calculate  new 


Algorithm  6  solves  for  refined  edge  flux.  It  accounts  for  all  within  cell 
scattering  through  inversion  of  the  scattering  source  as  was  done  for  SC.  The 
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beta  iteration  loop  defines  sonrce  distribntion  and  facilitates  calculation  of  cell 
transport  coefficients.  The  coupled  beta  edge  flux  iteration  loop  solves  for  within 
cell  scattering  and  accounts  source  distribution,  it  does  not  explicitly  couple  cells 
across  the  spatial  domain.  However,  if  beta  is  held  fixed  at  a  current  iteration 
estimate  edge  flux  can  be  calculated  from  equations  (89)  and  (90)  just  as  was 
done  for  SC.  These  equations  can  be  collapsed  with  flux  weights,  as  was  done  for 
SC,  using  equations  (197)  through  (202).  The  spatially  coupled  transport 
equations  are  then  readily  solved  as  in  chapter  2  using  equation  (172).  This 
suggests  the  following  iteration  scheme.  Algorithm  8  that  accounts  for  refined 
angle  within  cell  transport,  the  exponential  source  distribution,  and  particle 
transport  across  cells  in  the  spatial  domain. 
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!  obtain  an  initial  guess  for  ^  ■ 

inj 

!  obtain  an  initial  guess  for  [3^ 

!  calculate  cell  coupling  coefficients  for  ^  (^[3^ 


-(0) 

estimate  initial  flux  \|/  - 

calculate  cell  coupling  coefficients  for 

Km 


calculate  ^  f  ^ 


ealculate  K  f  ^ 


calculate 


-(/■+1)  {-{b) 

calculate  y/)^  ^  I 


calculate  S 


calculate  S 


calculate  new  P^  ^  1 5’^^- ,  5* Xj  j 
End  dos^i^D  i converged 


use  refined  angular  edge  flux  to  collapse  transport  coefficients 
calculate  y/ 

"H 

apportion  coarse  to  obtain  new  fine 

l—{b)  n  ^ 

End  Do  when  ,  yf  less  than  convergence  tolerance 


1™ _ ij-r,™  o.  A _ .,,1 _ _ 1  c _ _ i _ i  171, tj- _ j.: _ 


Algorithm  8  couples  the  angular  and  spatially  explicit  edge  flux  solutions 
in  a  scheme  that  accounts  for  infinite  scatters  from  cell  to  cell  and  within  a  cell 
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as  was  done  for  SC  but  adds  the  inner  loop  to  solve  for  y^-then  iterates  the  edge 
flux  and  beta  iteration  schemes  to  get  a  best  value  for  edge  flux  distribution. 

This  estimate  is  used  to  update  collapsed  transport  coefficients  which  are  used  to 
calculate  spatially  coupled  coarse  edge  flux  and  to  apportion  this  coarse  edge  flux 
back  into  fine  flux  components.  The  apportioned  edge  flux  updates  the  cell 
incoming  edge  flux  and  the  iteration  scheme  proceeds  again.  The  iteration 
scheme  of  Algorithm  8  effectively  accounts  for  infinite  particle  flights  in  a  single 
iteration.  Results  of  this  algorithm  implemented  as  Fortran  95  code  are 
discussed  in  the  next  chapter. 
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V.  EC  Experimental  Results 


After  showing  that  the  ASDI  method  met  my  evaluation  criteria:  robust  (reliably 
accurate  for  the  problems  intended),  efficient  (requires  few  iterations  to  converge), 
effective  (converges  rapidly  as  measured  by  compute  time)  with  SC,  I  implemented  the 
method  using  EC.  I  used  the  same  suite  of  experiments  to  test  the  methods  accuracy, 
effectiveness  and  computationally  efficiency  for  the  EC  experiments  as  I  did  for  SC. 

To  review,  the  first  experiment  had  an  optically  thick  homogeneous  material.  I 
varied  this  homogeneous  material  from  absorber  (low  scattering  ratio)  to  scatterer 
(scattering  ratio  nearly  one)  in  order  to  test  ASDI  against  SI  for  problems  which 
converged  readily  and  one  for  problems  which  converged  slowly. 

The  second  experiment  investigated  the  effect  of  periodic  material  discontinuity 
on  the  comparative  performance  of  ASDI  by  using  two  materials  of  the  same  dimension. 
There  was  an  emission  source  in  one  material  and  no  source  in  the  other  material.  Both 
materials  were  1  MFP  wide.  The  emitter  had  a  scattering  ratio  of  one.  The  non¬ 
emitter’s  scattering  ratio  varied  from  0.0  to  1.0  as  in  the  homogeneous  material.  This 
two  material  pattern  was  repeated  10  times  creating  a  periodic  discontinuity  in  the 
materials. 

I  again  investigated  the  effect  of  scattering  ratio,  spatial  refinement,  angular 
refinement  and  convergence  criterion  on  robustness,  effectiveness,  and  efficiency  for  both 
the  homogeneous  and  periodic  problems.  The  parameters  I  varied  in  turn  were: 
scattering  ratio,  spatial  refinement,  angular  refinement  and  convergence  tolerance.  I 
examined  these  parameters  for  a  homogeneous  then  a  two  material  problem.  Further,  I 
evaluated  ASDI  with  EC  over  a  complete  range  of  scattering  ratios  and  cross  sections, 
just  as  was  done  for  ASDI  with  SC  and  displayed  in  figures  26,  27,  and  28. 

As  I  did  with  SC,  I  tested  ASDTs  accuracy  by  checking  its  solution  against  a 
benchmark  to  ensure  it  met  the  convergence  tolerance  required.  The  benchmarks  were 
the  same  as  those  used  for  SC. 
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Problem  1  Optically  Thick  Homogeneous  Material 

The  first  test  problem  studied  is  a  100  cm  thick  homogeneous  slab  with  a  vacuum 
boundary  on  the  right,  a  symmetry  boundary  on  the  left,  a  uniform  isotropic  source, 
and  a  total  cross  section  of  1.0  cm'^ .  Figure  4,  although  shown  in  chapter  3  is 
redisplayed  in  this  chapter  for  the  reader’s  convenience. 


Figure  4  Problem  1 

The  first  experiment  tested  the  EC  ASDI  solution  robustness  as  scattering  ratio 
varied.  Figure  29  displays  two  plots.  These  are  the  symmetric  relative  difference  versus 
scattering  ratio  of  the  ASDI  solution  and  the  SI  solution  (ASDI_SI)  and  of  the  ASDI 
solution  and  the  SI  solution  after  being  given  the  ASDI  solution  (ASDI_SI-ASDI).  The 
convergence  tolerance  chosen  was  10  ®  for  both  SI  and  ASDI.  The  SI  Solution  after 
having  been  given  the  ASDI  solution  was  allowed  to  iterate  the  same  number  of  times 
that  it  took  to  arrive  at  its  converged  solution  without  the  ASDI  start  point.  For  this 
problem  this  is  between  30  and  500  iterations.  The  plot  shows  that  the  ASDI  solution 
never  deviates  from  this  solution  by  more  than  the  convergence  tolerance  required.  This 
demonstrates  ASDFs  reliable  accuracy  at  all  scattering  ratios.  At  a  scattering  of  0.6 
and  above  the  SI  solution  differs  by  more  than  the  chosen  convergence  tolerance 
demonstrating  Si’s  unreliable  solutions.  This  difference  increases  as  the  scattering  ratio 
increases.  These  plots  show  that  SI  produces  unreliable  answers  at  scattering  ratios  as 
low  as  0.6  but  ASDI  provides  reliably  accurate  answers  validated  by  Si’s  fixed  point 
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recognition.  In  fignre  6  (y45D/,  57  _y45'Z)/)  increases  from  a  low  valne  of  10  ^  and  a 

high  value  of  5.0  *  10”^ .  This  range  of  values  is  within  the  required  convergence 
tolerance. 


C  fraction 

Figure  29:  Symmetric  Relative  Difference  between  the  ASDI  solution  and  the  SI 
solution  as  scattering  ratio  varies.  Angular  quadrature  is  DE-8,  refinement  is  50 
(crAa;  =  2MFPs),  convergence  tolerance  is  10~®. 

Figure  30  shows  iteration  count  versus  scattering  ratio  (c)  with  an  angular 
refinement  of  8,  a  mesh  width  of  2  MFPs  and  a  convergence  tolerance  oflO^®.  As  in 
SC,  iteration  count  increases  with  scattering  ratio,  climbing  steeply  as  c  approaches 
unity.  This  again  demonstrates  the  classic  weakness  of  SI  for  problems  with  little  or  no 
absorption.  With  exponential  characteristic  spatial  quadrature  each  of  these  iterations 
requires  the  non-linear  root  finding  of  the  source  distribution  parameter  {Pni)  every 
direction  in  every  cell.  This  root  finding  can  be  computationally  expensive.  The  ASDI 
solution  converges  in  six  or  less  iterations  regardless  of  scattering  ratio,  demonstrating 
its  effectiveness  even  for  problems  with  little  or  no  absorption. 
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Figure  30  Plot  of  iteration  count  versus  scattering  ratio.  Angular  quadrature  is  DE-8, 
refinement  is  50  (Ax  =  2MFPs  ),  convergence  tolerance  is  10~®. 

As  pointed  out  previously  the  low  iteration  count  of  the  ASDI  method  with  EC 
displayed  in  Figure  30  is  not  a  practical  improvement  over  SI  unless  computational  cost 
is  also  reduced.  Figure  31  displays  compute  time  versus  scattering  ratio.  The  amount 
of  compute  time  needed  by  the  ASDI  method  is  less  than  0.3  seconds  regardless  of  the 
scattering  ratio.  Even  when  SI  converges  more  rapidly  than  ASDI  and  provides  reliably 
accurate  solutions,  i.e.  at  scattering  ratios  below  0.6,  ASDI  compute  time  is  the  same 
order  of  magnitude  as  SI.  The  plot  confirms  that  the  beta  feedback  mechanisms  of  the 
ASDI  method  with  EC  do  not  converge  slowly  as  the  scattering  ratio  goes  to  one.  This 
demonstrates  that  the  method  does  not  require  excessive  compute  time  with  the  EC 
spatial  quadrature  at  least  for  the  homogeneous  material. 
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Figure  31  Plot  of  compute  time  (seconds)  versus  scattering  ratio.  Angular 


quadrature  is  DE-8,  refinement  is  50  {a Ax  =  2MFPs),  convergence  tolerance  is  10 


The  results  shown  in  Fignre  30  were  obtained  with  a  fixed  refinement  of  50  mesh  cells 
{Ax  =  2.0cm),  convergence  tolerance  fixed  at  10^®,  and  an  angnlar  qnadratnre  fixed  at 

n  =  8 .  Spatial  refinement,  angnlar  refinement,  and  convergence  tolerance  are  also 
parameters  that  are  expected  to  have  an  impact  on  robnstness,  effectiveness,  and 
efficiency  of  the  ASDI  method  nsed  with  the  EC  spatial  qnadratnre.  The  next  series  of 
plots  investigates  the  effect  changing  these  parameters  have  on  ASDI  performance  with 
the  EC  spatial  qnadratnre.  These  three  parameters  were  stndied  with  a  scattering  ratio 
of  0.98  making  the  homogeneons  material  a  good  scatterer. 

Fignre  32  shows  the  relative  difference  between  SI,  ASDI  and  an  estimated  best 
solntion  which  is  SI  after  being  given  the  ASDI  solntion  and  allowed  to  iterate  enongh 
times  to  confirm  a  fixed  point.  The  difference  between  the  ASDI  solntion  and  the  SI- 
ASDI  solution  is  less  than  the  convergence  tolerance  for  all  refinements  whereas  the 
difference  between  SI  and  the  SI-ASDI  answers  is  greater  than  the  convergence 
tolerance.  This  is  because  SI  converges  falsely  in  this  optically  thick  problem  with  little 
or  no  absorption. 
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Figure  32  Plot  of  Symmetric  Relative  Difference  £srd  between  ASDI  solution  and  SI 
solutions  for  as  cell  mesh  is  refined.  Angular  quadrature  is  DE-8,  scattering  ratio  is  1.0, 

convergence  tolerance  is  10  ® . 


SI  recognition  of  the  ASDI  solntion  as  a  fixed  point  for  all  spatial  refinements 
provides  convincing  evidence  of  the  accnracy  of  its  solntion.  I  was  able  to  farther  test 
this  accnracy  becanse  an  analytic  solntion  was  available  for  problem  1.  I  compnted  and 
compared  the  ASDI  convergence  rate  with  the  known  convergence  rate  of  EC.  This  was 
done  for  a  DE-4  angnlar  qnadratnre  and  a  scattering  ratio  of  1.0  with  spatial  refinement 
varying  from  10  (cell  width=10cm)  throngh  1000  (cell  width  =  1.25xl0~^cm).  I  did 

not  expect  the  ASDI  solntions  be  the  same  as  the  analytic  solntion  for  coarse  spatial 
meshes  becanse  the  spatial  qnadratnre  is  a  nnmerical  approximation.  However,  becanse 
EC  is  known  to  be  fonrth  order  convergent  in  space,  the  method  can  be  checked  against 
the  analytic  solntion  to  determine  if  the  order  of  spatial  convergence  is  fonr.  Fignre  33 
shows  symmetric  relative  difference  (  £srd  )  between  the  ASDI  solntion  and  an  exact 
solntion  plotted  against  the  spatial  mesh  refinement  factor  (Rp)  for  a  scattering  ratio  of 
1.0.  SsRD  appears  as  a  straight  line  on  both  plots  whose  slope  is  3.86.  A  line  with  a 
slope  of  fonr,  the  convergence  order  of  EC,  has  been  overlaid  on  the  data  plot.  This 
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agreement  of  ASDI’s  convergence  rate  with  the  known  EC  convergence  rate  for  a 
scattering  ratio  of  1.0  combined  with  SI  recognition  of  ASDI  solutions  as  fixed  points 
demonstrate  the  method’s  robustness.  The  preservation  of  EC  fourth  order  convergence 
with  spatial  refinement  provides  strong  evidence  that  the  ASDI  method  is  not  changing 
the  spatial  convergence  properties  of  the  EC  spatial  quadrature. 


Figure  33:  Plot  of  Symmetric  Relative  Difference  between  the  ASDI  solntion  and 
an  analytic  solntion  as  cell  mesh  is  refined  between  10  MFPs  and  0.1  MFPs.  Angular 
quadrature  is  DE-8,  scattering  ratio  is  1.0,  convergence  tolerance  is  10~® . 


Figure  34  shows  the  iteration  count  of  the  ASDI  and  SI  methods  as  the  mesh  is 
refined.  SI  iteration  remains  flat,  at  500  iterations,  ASDI  iteration  count  is  seven  or  less 
iterations  through  mesh  refinements  of  64.  This  indicates  that  the  number  of  iterations 
required  for  convergence  is  not  dependent  on  the  mesh  refinement  for  ASDI  even  when 
applied  to  the  non-linear  EC  spatial  quadrature. 
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Refinement  Factor 

Figure  34  Plot  of  iteration  count  as  cell  mesh  is  refined.  Angular  quadrature  is  DE-8, 
scattering  ratio  is  1.0,  convergence  tolerance  is  . 


Figure  35  shows  the  compute  times  required  for  the  SI  and  ASDI  solutions  as  the 
spatial  mesh  is  refined.  The  spatial  refinement  increases  compute  time  for  both 
methods.  ASDI  requires  less  compute  time  than  SI  through  a  refinement  of  32. 
Although  the  ASDI  method  requires  more  compute  time  than  SI  for  refinements  above 
32.  This  behavior  was  also  observed  for  ASDI  applied  to  SC  in  figure  12.  The  plots 
indicate  that  for  a  fine  enough  mesh  source  iteration  will  converge  faster  than  ASDI. 
However,  SI  will  not  converge  to  reliable  and  accurate  solutions.  The  additional  time 
used  to  compute  ASDI  approximations  results  accurate  solutions.  The  plot  indicates 
that  ASDI  will  not  require  disproportionately  more  time  than  SI  for  refined  spatial 
meshes  and  will  converge  to  the  correct  solution.  The  results  observed  for  EC  are 
similar  to  those  obtained  for  SC  indicating  the  beta  subroutine  is  working  as  intended. 
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Figure  35  Plot  of  compute  time  required  for  the  ASDI  and  SI  solutions  as  cell  mesh  is 
refined.  Angular  quadrature  is  DE-8,  scattering  ratio  is  1.0,  convergence  tolerance  is 

10"®. 


The  convergence  tolerance  nsed  to  test  previons  problems  was  10"®  on  angnlar 
flnx.  This  is  a  fairly  tight  convergence  tolerance  and  it  is  snitable  for  most  engineering 
applications.  However,  some  applications  may  reqnire  tighter  tolerances.  The  effect  of 
tightening  convergence  tolerance  from  10"® to  10"^  on  symmetric  relative  difference, 
iteration  connt,  and  compnte  time  is  shown  next.  Fignre  36  shows  that  the  accnracy  of 
the  converged  SI  solntion  does  not  meet  the  accnracy  reqnired  by  the  specified 
convergence  tolerance  regardless  of  how  tight  that  tolerance  is,  whereas  the  ASDI 
method  continnes  to  provide  reliably  accnrate  solntions  for  any  tolerance  withont  ronnd 
off  error  or  instability.  Examination  of  Fignre  36  and  Fignre  37  reveals  that  this 
accnracy  is  achieved  with  modest  increase  in  iterations  reqnired  and  compnte  time  for 
the  ASDI  method  applied  to  EC.  These  results  are  similar  to  those  obtained  for  SC  and 
displayed  in  Figures  13,  14,  and  15. 
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Convergence  Tolerance 

Figure  36  Plot  of  symmetric  relative  difference  versus  convergence  tolerance.  Angular 
quadrature  is  DE-8,  scattering  ratio  is  1.0,  cell  size  is  2  MFPs. 


Convergence  T  olerance 


Figure  37  Plot  of  iterations  versus  convergence  tolerance  for  ASDI  and  SI.  Angular 
quadrature  is  DE-8,  scattering  ratio  is  1.0,  cell  size  is  2  MFPs. 
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Figure  38  Plot  of  compute  time  versus  convergence  tolerance  for  ASDI  and  SI. 
Angular  quadrature  is  DE-8,  scattering  ratio  is  1.0,  cell  size  is  2  MFPs. 


The  effect  of  increasing  angnlar  refinement  for  ASDI  applied  to  EC  is  shown  in 
Fignre  39  throngh  Fignre  41.  These  plots  demonstrate  that  ASDI  continnes  to  provide 
reliably  accnrate  solntions  withont  significantly  increasing  iteration  connt  and  with 
practical  compnte  times  as  the  nnmber  of  ordinates  in  the  angnlar  qnadratnre  is 
increased.  SI,  as  expected,  remains  inaccnrate  and  increases  compntational  cost  as  the 
angnlar  qnadratnre  is  refined.  As  with  ASDI  applied  to  SC,  ASDI  applied  to  EC 
converges  with  the  least  iteration  and  most  rapid  compnte  times  when  the  angnlar 
qnadratnre  has  only  two  directions.  ASDI  achieves  a  qnick  solntion  for  this  coarse 
angnlar  qnadratnre  becanse  the  global  spatial  rontine  solves  the  same  problem  on  the 
first  iteration  as  the  fine  angle  rontine.  After  the  first  iteration  both  edge  flnx  and  beta 
are  nearly  correct  and  converge  rapidly  with  snccessive  iteration.  Unfortnnately  such  a 
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coarse  angular  quadrature  is  not  likely  to  meet  most  engineering  needs.  Figure  39  and 


Figure  40  display  the  same  convergence  reliability  characteristics  and  iteration  count 
versus  number  of  directions  for  EC  as  was  observed  with  SC  in  Figure  16  and  Figure  17. 
This  is  encouraging  because  the  increased  angular  refinement  requires  increasing  the 
number  of  calls  to  the  beta  root  finding  algorithm.  The  increased  number  of  beta  root 
finding  problems  does  not  affect  reliability  or  significantly  increase  compute  times. 
However,  comparison  of  Figure  41,  which  shows  ASDI-EC  compute  time  versus  angular 
refinement,  and  Figure  18,  which  shows  ASDI-SC,  compute  time  versus  angular 
refinement,  demonstrates  that  compute  time  increases  more  steeply  for  ASDI-EC  than  it 
does  for  ASDI-SC.  The  steeper  increase  in  compute  time  for  EC  occurs  because  beta 
must  be  found  for  every  direction  is  every  cell.. 


Figure  39  Plot  of  symmetric  relative  difference  versus  angular  quadrature. 
Convergence  tolerance  is  10  ®  scattering  ratio  is  1.0,  cell  size  is  2  MFPs. 
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Figure  40  Plot  of  iterations  versus  angular  quadrature  for  ASDI  and  SI.  Convergence 
tolerance  is  10~®  scattering  ratio  is  1.0,  cell  size  is  2  MFPs. 


Figure  41  Plot  of  compute  time  versus  angular  quadrature  for  ASDI  and  SI. 
Convergence  tolerance  is  10  ®  scattering  ratio  is  1.0,  cell  size  is  2  MFPs. 
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This  section  demonstrated  the  strength  of  the  ASDI  method  with  the  EC  spatial 
qnadratnre  for  a  homogeneons  material.  The  most  important  characteristic  examined 
was  accnracy  for  each  of  the  parameters  stndied.  Each  experiment  showed  that  ASDI 
applied  to  EC  was  reliably  accnrate  within  the  convergence  tolerance  reqnired  whereas 
SI  was  not.  The  application  of  ASDI  to  the  non-linear  EC  method  demonstrated  very 
similar  resnlts  to  the  method’s  application  to  SC.  This  confirmed  the  hypothesis  that 
beta  conld  be  fonnd  and  fixed  for  each  iteration  resnlting  in  an  algorithm  that  was 
similar  to  that  which  obtained  ontstanding  result  with  SC.  The  EC  algorithm 
essentially  adds  a  root  finding  subroutine.  This  result  shows  that  the  ASDI  method  is 
not  limited  to  a  linear  spatial  quadrature.  It  indicates  a  more  general  application.  The 
method’s  reliable  accuracy,  low  iteration  count,  and  fast  compute  times  across  all 
parameters  studied  demonstrate  that  the  method  is  robust,  effective  and  efficient  with 
EC  for  this  homogeneous  problem.  The  next  section  examines  the  method’s 
performance  when  applied  to  a  heterogeneous  material. 

Periodic:  Two  Regions  Repeated  10  Times 

The  second  test  problem  investigated  is  the  same  problem  investigated  second  with  SC. 
It  is  a  slab  with  a  symmetry  boundary  on  the  left  side  and  a  vacuum  boundary  on  the 
right  side.  Two  materials  of  I  cm  cell  width  are  placed  side  by  side.  This  two  material 
pattern  is  repeated  10  times  for  a  total  length  of  20  cm.  Material  I  had  a  total  cross 
section  of  1.0  cm"  ,  a  scattering  ratio  of  1.0,  and  a  uniform  source  of  1.0  cm"  .  Its 
material  properties  remain  fixed  for  all  experiments.  Material  II  had  a  total  cross 
section  of  1.0  cm'^ ,  a  baseline  scattering  ratio  of  1.0,  and  no  source.  Material  II’s 
parameters  were  varied  during  the  experiments.  A  diagram  of  these  two  materials  was 
first  shown  in  chapter  three  with  Figure  19.  It  is  redisplayed  here  for  the  reader’s 
convenience. 
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Figure  19  Problem  2 

As  with  SC  the  first  set  of  experiments  for  problem  two  tested  the  ASDI 
method’s  accuracy,  effectiveness  and  efficiency  versus  scattering  ratio  fraction  in 
Material  11.  Scattering  ratio  fraction  was  varied  from  0.0  to  1.0.  Angular  refinement  for 
this  experiment  was  8,  refinement  factor  was  1  and  convergence  tolerance  was  10 

Accuracy  is  the  first  parameter  presented.  Figure  42  displays  the  symmetric 
relative  difference  between  ASDI  and  SI  (ASDI,SI_ASDl))  and  between  ASDI  and 
SI  given  the  ASDI  solution  (ASDI,SI_ASDl)).  The  plots  show  that  ASDI  applied 
to  EC  provides  reliably  accurate  solutions  for  all  scattering  ratios  while  SI  provides 
reliably  accurate  solutions  only  through  scattering  ratios  of  0.6.  The  reliably  accuracy 
of  Figure  42  for  EC  is  similar  to  the  reliable  accuracy  of  Figure  20  for  SC.  This  shows 
that  the  beta  convergence  does  not  effect  reliable  accuracy  of  the  method,  even  in 
problems  that  are  not  homogeneous. 
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Periodic  EC 


Figure  42:  Symmetric  Relative  Difference  egjuj  between  the  ASDI  solntion  and  the  SI 
solntion  as  scattering  ratios  vary.  Angnlar  qnadratnre  is  DE-8,  convergence  tolerance 

is  10  ® ,  refinement  is  50. 


The  effectiveness  of  the  method  was  tested  next.  Figure  43  shows  the  method 
iteration  count  versus  scattering  ratio  fraction  for  material  11.  This  plot  shows  that 
ASDI  converges  in  less  than  six  iterations  regardless  of  scattering  ratio  demonstrating 
its  effectiveness,  even  when  applied  to  EC,  for  heterogeneous  materials  with  little  or  no 
absorption.  These  results  are  similar  to  those  obtained  from  SC  and  demonstrate  that 
beta  convergence  did  not  degrade  effectiveness. 
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C  fraction 

Figure  43:  Plot  of  iteration  count  versus  scattering  ratio  fraction.  Angular  quadrature 
is  DE-8,  convergence  tolerance  is  10  ® ,  refinement  is  50. 

Efficiency  of  the  method  is  tested  next.  Fignre  44  shows  compnte  time  versns 
scattering  ratio  fraction  for  material  11.  This  fignre  shows  that  ASDI  compnte  time  is 
less  than  0.1  seconds  for  all  scattering  ratios  tested.  Comparison  of  the  ASDI  and  SI 
compnte  times  demonstrate  that  ASDI  is  more  compntationally  efficient  than  SI  even 
when  material  II  is  a  strong  absorber  and  SI  converges  rapidly.  The  fignre  also  shows 
that  as  scattering  ratio  approaches  one  in  material  II  SI  compnte  time  climbs  steeply 
bnt  ASDI  compnte  times  do  not.  These  results  are  similar  to  the  compute  times 
displayed  in  Figure  22  for  SC.  The  results  demonstrate  that,  although  the  beta 
convergence  loop  was  added  to  the  ASDI  algorithm  for  EC,  the  algorithm  does  not 
require  excessive  compute  time. 


5-17 


Periodic  EC 


Figure  44:  Plot  of  compute  time  (seconds)  versus  scattering  ratio.  Angular  quadrature 
is  DE-8,  convergence  tolerance  is  10  ® ,  refinement  is  50. 


The  next  set  of  experiments  for  problem  two  demonstrate  method  accnracy, 
effectiveness  and  efficiency  while  refinement  factor  varies  for  material  11.  Scattering 
ratio  for  material  II  was  fixed  at  1.0.  This  tested  ASDI  in  materials  with  little  or  no 
absorption.  A  material  discontinnity  exists  even  thongh  cross  sections  and  scattering 
ratios  are  the  same  becanse  intrinsic  sonrces  are  emitted  only  in  material  I.  The 
refinement  factor  for  material  II  was  varied  between  1  and  66. 

Accnracy  was  tested  first.  Fignre  45  displays  the  symmetric  relative  difference 
between  ASDI  and  SI  (ASDI,SI_ASDl))  and  between  ASDI  and  SI  given  the  ASDI 
solntion  (£'§j^j)(ASDI,SI_ASDI)).  The  plots  show  that  ASDI  provides  reliably  accnrate 
solntions  for  all  refinement  factors  bnt  SI  does  not. 
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Figure  45:  Plot  of  Symmetric  Relative  Difference  between  ASDI  solution  and  SI 
solutions  as  cell  mesh  is  refined.  Angular  quadrature  is  DE-8,  scattering  ratios  are 
baseline  values,  convergence  tolerance  is  10~® . 
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Effectiveness  was  next  tested.  Fignre  46  shows  that  ASDI  converges  on  a  solntion  in 
less  than  8  iterations  for  a  refinement  factor  of  1  growing  slightly  as  refinement  factor  is 
increased  then  remaining  constant  for  refinement  factors  above  16.  SI  requires  many 
more  iterations  than  ASDI  and  provides  unreliable  solutions  for  these  refinement 
factors. 
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Figure  46  Plot  of  iteration  count  as  cell  mesh  is  refined.  Angular  quadrature  is  DE-8, 
scattering  ratios  are  baseline  values,  convergence  tolerance  is  10~® . 


Computational  efficiency  was  next  tested.  Figure  47  demonstrates  that  ASDI 
requires  less  compute  time  than  SI  regardless  of  mesh  size  for  the  baseline  scattering 
ratios  selected.  However,  as  spatial  mesh  is  refined  the  compute  time  of  ASDI  is  nearly 
the  same  as  the  compute  time  of  SI.  The  relatively  comparable  compute  times  of  SI 
and  ASDI  do  not  indicate  that  SI  is  just  as  efficient  as  ASDI  for  these  fine  spatial 
meshes  because  SI  does  not  provide  reliably  accurate  answers.  ASDI  does  not  require 
excessive  compute  time  regardless  of  refinement  factors. 
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Figure  47: 


Plot  of  compute  time  as  cell  mesh  is  refined.  Angular  quadrature  is  DE-8, 


scattering  ratios  are  baseline  values,  convergence  tolerance  is  10 
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Accuracy,  effectiveness  and  compntational  efficiency  were  fnrther  tested  with 
fixed  scattering  ratio  and  refinement  factor  for  varying  convergence  tolerance  and 
angnlar  refinement  with  similar  results.  As  in  problem  1  these  results  indicate  that 
ASDI  remains  reliably  accurate,  effective  and  efficient  as  these  parameters  vary.  This  is 
the  same  conclusion  drawn  for  SC  and  provides  strong  confirmation  that  the  non¬ 
linearity  of  the  EC  spatial  quadrature  algorithm  does  not  prevent  applying  the  ASDI 
algorithm  to  EC  which  is  what  I  set  out  to  do. 

These  experiments  demonstrate  the  strength  of  the  ASDI  method  for  a  two 
region  periodic  material.  The  experiments  show  that  ASDI  was  always  reliably  accurate 
and  that  the  parameters  with  the  greatest  impact  on  iteration  count  and  compute  time 
were  scattering  ratio  and  refinement  factor.  The  experiments  presented  kept  cross 
sections  constant  but  varied  spatial  refinement.  In  order  to  fully  stress  material 
discontinuity  in  this  periodic  problem  I  next  examined  the  impact  of  varying  both 
scattering  ratio  and  cross  section  simultaneously  in  material  II,  just  as  I  did  with  SC.  I 

_ n  _ 1 

varied  scattering  ratio  fraction  from  0.0  to  1.0  and  cross  section  from  10  cm  to 
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10  cm  ^ .  This  range  of  cross  sections  is  not  the  same  as  the  range  of  cross  sections  nsed 
for  SC  (i.e.  10~^cm”^  to  10^  cm~^).  That  is  becanse  the  exponential  characteristic 
method  as  cnrrently  implemented  can  be  nsed  for  optical  thicknesses  of  abont  80. 

The  limitation  in  cross  section  resnlts  from  nnmerically  poor  conditioning  in  the 
calculation  of  EC  transport  coefficients.  It  is  not  a  limitation  of  the  ASDI  method.  I 
examined  the  impact  of  varying  these  two  parameters  on  accuracy,  effectiveness  and 
efficiency  at  a  refinement  factor  of  1. 

A  3D  plot  of  symmetric  relative  difference  between  ASDI  and  SI 
( ^SRD  (ASDI, si)  )  and  between  ASDI  and  SI  given  the  ASDI  solution 
( ^'SRD  (ASDI,SI_ASDI)  )  is  shown  in  Figure  48.  This  plot  shows  that  ASDI  is  reliably 
accurate  across  the  range  of  material  properties  tested.  The  plot  shows  a  spike  in 
SsRD  {ASDI,SI  ASDI^dX  cross  sections  of  approximately  1 0 .  Since  this  spike  is  in 
the  region  that  is  below  the  convergence  tolerance  that  was  required  (10~^)  I  draw  no 
conclusions.  The  spike  indicates  that  cross  section  values  greater  than  lOcm”^  might  be 
unreliable. 
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Figure  48:  3D  Plot  of  ASDI  symmetric  relative  difference  between  ASDI  and  SI 
( ^SRD  ( ASDl,Sl_ASDl) )  and  between  ASDI  and  SI  given  the  ASDI  solution 
(f’sRD  (ASD1,S1_ASD1)  )  as  scattering  ratio  and  cross  section  Vciry  .  Angular  quadrature 
is  DE-8,  convergence  tolerance  is  ,  refinement  factor  is  1. 

A  3D  plot  of  ASDI  iteration  count  versus  cross  section  and  scattering  ratio  is 
shown  in  Figure  49. 
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Figure  49:  3D  Plot  of  ASDI  iteration  count  as  scattering  ratio  and  cross  section  vary 


Angular  quadrature  is  DE-8,  convergence  tolerance  is  10  ,  refinement  factor  is  1. 


This  figure  demonstrates  that  ASDI  is  not  sensitive  to  scattering  ratio  for  cross  sections 
less  than  \cm~^ .  For  cross  sections  less  than  IcnT^  ASDI  takes  7  iterations  or  less  to 
converge.  Iteration  connt  increases  sharply  near  scattering  ratios  of  one  and  cross 
sections  of  lOcm”^.  Fnrther  investigation  is  needed  to  determine  the  reason  for  this 
increase  in  iteration  count.  It  was  not  present  with  ASDI  applied  to  SC  as  displayed  in 
Figure  27. 

As  with  SC  a  low  iteration  count  is  not  useful  if  computational  cost  is 
prohibitive.  A  comparison  of  the  time  required  for  ASDI  and  SI  to  solve  the  problem  is 
shown  in  Figure  50.  The  ratio  of  the  ASDI  compute  time  to  SI  compute  time  is  shown 
in  the  figure,  if  this  ratio  is  less  than  1.0  then  ASDI  takes  less  time  than  SI.  This  is  the 
case  for  cross  sections  less  than  10  cm  .  The  figure  shows  that  SI  takes  one  tenth  the 
compute  time  of  ASDI  for  cross  sections  greater  than  10  cm  with  scattering  ratios  of 
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0.1.  In  this  region  SI  converges  in  9/100  second  and  ASDI  converges  in  60/100  second. 

SI  converges  very  rapidly  becanse  the  problem  is  absorptive  in  more  than  half  problem 

material.  However,  SI  converges  falsely.  This  false  convergence  occnrs  becanse  half  the 

problem  (material  one)  has  little  or  no  absorption  even  when  the  scattering  ratio  of 

material  two  is  nearly  zero.  In  general,  SI  takes  less  compnte  time  than  ASDI  when 

scattering  ratios  are  less  than  0.6  and  cross  sections  are  greater  than  10  cm  .In  this 

domain  SI  converges  in  nearly  no  time  at  all.  Even  in  these  problems  ASDI  compute 

times  are  practical  and  competitive  with  SI  compute  times.  More  importantly,  ASDI 

provides  reliably  accurate  answers  across  the  full  range  of  cross  section  and  scattering 

ratio  parameters  just  as  it  did  with  SC.  Further  research  is  required  to  determine  if  the 

—2  —1 

greater  compute  times  required  for  cross  sections  greater  than  10  cm  result  from  the 
calculation  of  EC  transport  coefficients. 
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Figure  50:  3D  Plot  of  the  ratio  ASDI  compute  time  to  SI  Compute  time  as  scattering 
ratio  and  cross  section  vary  .  Angular  quadrature  is  DE-8,  convergence  tolerance  is 

10~® ,  refinement  factor  is  1. 

These  experiments  demonstrate  that  ASDI  is  reliably  accnrate  for  all  the 
parameters  stndied  and  therefore  robnst.  The  method  dramatically  rednces  iteration 
count  and  required  compute  time  for  diffusive  problems.  The  experimental  results  show 
that  ASDI  is  a  useful  and  practical  iteration  scheme  for  any  material  properties.  The 
method  now  allows  researchers  to  examine  the  computation  of  EC  transport  coefficients 
in  problems  that  are  optically  thick  with  little  or  no  absorption. 
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VI.  Summary  and  Conclusions 


Accurate,  and  robust  positive  spatial  quadrature  schemes  used  in  discrete  ordinates 
methods,  such  as  EC,  provide  physically  meaningful,  non-negative  fluxes  given  non¬ 
negative  incident  fluxes,  non- negative  emission  sources,  and  non- negative  scattering 
cross  sections.  In  1969  K.D.  Lathrop  eloquently  stated  why  positive  spatial 
quadratures  are  needed.  He  said  that  in  addition  to  numerical  difficulties  there  are 
“psychological  problems”  with  negative  fluxes.  The  user  who  understands  the  transport 
equation  and  not  just  the  numerical  solution  procedure  knows  that  there  is  no  such 
thing  as  a  negative  angle  integrated  flux,  and  rapidly  becomes  cynical  about  the 
effectiveness  of  a  program  which  produces  negative  numbers  {12  :476)”.  Mathews  et  all 
have  been  developing  the  EC  spatial  quadrature  since  the  early  1990’s  which  is  positive 
and  approaches  fourth  order  in  its  spatial  convergence.  To  date  this  method  has  been 
difficult  to  implement  in  problems  where  scattering  ratios  were  nearly  one  (7:36, 

11:165).  EC  and  other  more  conventional  spatial  quadratures  based  on  source  iteration 
are  impractical  for  highly  diffusive,  optically  thick  problems.  The  objective  of  this 
effort  was  to  develop  an  accurate  and  efficient  scheme  to  rapidly  converge  EC.  This  has 
been  done  for  slab  geometry.  This  contribution  makes  it  possible  to  use  the  exponential 
characteristic  method,  and  similar  methods,  for  real  materials,  even  for  notoriously  slow 
transport  computations  involving  one  or  a  few  groups  describing  the  thermal  neutron 
energy  range  from  about  0  to  1  eV.  (2:83).  I  originally  explored  synthetic  acceleration 
as  the  method  to  rapidly  converge  these  quadratures.  Although  successful  for  most 
problems  this  technique  diverged  for  problems  with  sharp  material  discontinuities.  By 
taking  advantage  of  the  lessons  learned  from  synthetic  acceleration,  I  found  it  possible 
to  formulate  a  new  transport  method  that  more  directly  solves  for  fine  angular 
resolution  flux  and  spatially  coupled  coarse  angular  resolution  flux.  The  fine  angular 
quadrature  accounts  for  the  contribution  of  infinite  number  of  particle  flights  to 
scattered  source  within  a  single  cell.  The  coarse  angular  quadrature  accounts  for  the 
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contribution  of  an  infinite  number  of  particle  flights  to  scattered  source  across  the 
spatial  domain. 

The  introduction  of  two  transport  techniques,  one  producing  full  angle  coupling 
for  a  fine  quadrature  and  one  producing  full  spatial  coupling  without  iteration  are  both 
major  contributions  of  this  research.  The  ability  to  couple  these  methods  and  iterate  on 
transport  coefficients  for  both  EC  and  SC  instead  of  source  iteration  removes  the  need 
to  accelerate  optically  diffuse  problems  per  se.  The  transport  method  provides  the 
estimate  of  infinite  particle  flights.  The  ASDI  method  provides  an  alternative  to 
conventional  source  iteration  at  the  expense  of  a  larger  linear  algebra  problem.  It  has 
the  demonstrated  advantage  robustness  (reliably  accurate  for  the  problems  it  is  designed 
for)  effective  (requires  few  iterations)  and  efficient  (requires  practical  compute  times). 

The  coupling  of  the  two  transport  methods  worked  surprisingly  well.  The 
discovery  that  EC  converges  rapidly  with  this  method  and  can  now  be  applied  to 
diffusive  problems  should  lead  to  renewed  interest  in  the  EC  spatial  quadrature  within 
the  transport  community. 

The  use  of  flux  weights  to  project  between  coarse  and  fine  angular  quadratures  is 
general  with  respect  to  the  spatial  quadrature  chosen.  Extension  to  other  positive 
spatial  quadratures  is  immediate.  This  should  further  generate  interest  within  the 
community  for  the  method  developed.  Spatial  quadratures  that  do  not  preserve  positive 
flux  require  further  research  regarding  flux  weighting.  I  did  not  derive  a  method  to 
calculate  flux  weights  when  faced  with  negative  currents.  Flux  weights  as  employed  in 
the  method  approximate  flux  distributions  at  cell  edges.  These  flux  weights  conserve 
angular  information  while  collapsing  to  a  coarse  quadrature  that  readily  captures 
diffusive  behavior.  The  concept  of  flux  weights,  and  cell  edge  flux  distribution  should 
find  application  in  a  wider  set  of  transport  applications. 

The  method  has  been  implemented  in  FORTRAN-95  on  a  PC,  demonstrating  the 
computational  practicality  of  the  approach  in  slab  geometry  with  isotropic  cross 
sections.  The  method  as  derived  is  immediately  extendable  to  anisotropic  scattering 
cross  sections.  Its  extension  to  multiple  dimensions  will  entail  flux  weighting  on  four 
cell  edges  vice  two  and  an  efficient  solver  for  the  spatially  coupled  coarse  angular 
quadrature.  The  minimum  bandwidth  of  the  coefficient  matrix  which  readily  inverts  in 
ID  will  require  research  to  develop  an  efficient  2D  solver.  However,  the  demonstrated 
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ability  to  cast  a  fine  angle  problem  into  an  eqnivalent  coarse  angle  problem  should  make 
the  effort  in  2D  tractable. 

The  fine  angle  routine  algorithm  6  is  readily  parallelizable.  The  routine  uses 
iteration  edge  flux  as  an  estimate  for  each  spatial  cell.  Substantial  improvement  in 
performance  can  be  obtained  by  parallelizing  this  subroutine.  Further,  profiling  reveals 
that  60%  of  computational  effort  (time)  is  spent  with  this  routine.  Hence  significant 
reduction  in  the  ASDI  method’s  overall  computational  efficiency  would  result  from  this 
parallelization.  The  method  as  developed  lends  itself  to  adoptive  spatial  and  angular 
meshing.  There  is  good  reason  to  believe  that  a  problem,  such  as  the  periodic  horizontal 
interface  tested,  could  use  a  coarse  angular  and  spatial  mesh  in  regions  without  material 
discontinuity  and  a  finer  mesh  in  regions  that  needed  to  capture  angular  or  steaming 
behavior.  Research  in  this  area  is  straightforward. 

The  code  written  for  this  research  is  not  a  production  code.  A  graphical  user 
interface  to  obtain  user  input  defining  a  flexible  range  of  problem  parameters  is  needed. 
This  user  interface  should  also  interface  with  available  cross  section  computational  tools 
such  as  the  NJOY  suite  of  algorithms  and  Gerts’  PAX  cross  sections.  Nevertheless,  the 
practicality  of  my  method  has  been  demonstrated  successfully  for  EC  which  is  what  I 
set  out  to  do. 
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Appendix  A 

This  appendix  demonstrates  that  j)  meets  the  reqnirements  for  a 

distance  fnnction  and  that  it  can  readily  be  applied  to  vectors  in  FORTRAN  90. 

Al.  x)as  a  distance  function 

A  non-empty  ,X,  set  together  with  a  ‘distance  fnnction’  i^^y)  is 
said  to  form  a  metric  space  provided  that: 


^5i;zj(T^)  =  0  iff 

x  =  y 

(1) 

^SRD  (b?  ^  ^SRD  (^’  b)  —  0 

Vx,y  eX 

(2) 

(y ,  x)  =  (x,  z)  +  Sss,o  {y,  z) 

yx,y,z^X 

(3) 

The  distance  fnnction  ^^^^(y,x)is 


^SRD{x^y)  =  ^  (4) 

The  nnmerator  in  eqnation  (4) 

2|x-y|  =  0  (5) 

is  zero  if  and  only  if  and  only  if 

|x-y|  =  0  (6) 
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This  occurs  only  if  x  =  y.  The  denominator  in  equation  (4)  is  only  zero  if  x  and 
y  are  zero  in  which  case  the  distance  between  them  is  zero.  Equation  (4)  is  then 
zero  only  if  equation  (6)  is  zero.  This  occurs  if  only  when  x=y  meeting 
requirement  one. 

If  X  is  not  equal  to  y,  and  x  >0,  and  y>0  then 


and 


Imax 


min 


=  2, 


(Ixl+lyl)  (kl+lrl) 

''I  I  I-'  I/max  ''I  I  1-^  I/max 


which  meets  requirement  two. 


If  0  <  X  <  z  <  y  then 


F-r|  _ 

|jc|+|z|  z+x 


and 


ly-^l  _ 

|y|+|z|  z+y  ■ 


Adding  equation  (10)  to  equation  (9)  results  in 


(7) 


(8) 


(9) 


(10) 
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(11) 


2-x  j_ 

z+x  z+y  x+y  * 

Multiplying  each  of  the  terms  by  the  appropriate  denominator  results  in 

(z  -x){z  +  y){x  +  y)  +  {y-  z){z  +  x){x  +  y)>{y-x){z  +  x){z  +  y) .  (12) 

Subtracting  the  right  side  from  the  left  side  of  equation  (12)  and  simplifying 
results  in 

{x-y){x-z){y-z)>0  (13) 

which  is  always  true  based  on  the  initial  condition  0  <  x  <  z  <  y  .  The  same  result 
can  be  shown  if  x,  y  or  z  are  swapped.  This  demonstrates  that  the  distance 
function  meets  the  condition  three. 
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A.2.  Application  of  j)  to  Vectors 

The  distance  function  j)  defined  by  equation  (4)  can  be  applied  to 

vectors  x,y  resulting  in  a  vector  of  distances  Ssrd  [x,y^  .  An  element  of  this 
distance  vector  is 

^SRDi  b)  “  ^SRD  ( T  ’  b,'  )  ■  (14) 


In  my  algorithm  I  am  interested  in  ensuring  that  no  element  of  the 
distance  vector  ssrd  is  greater  than  a  convergence  tolerance.  Therefore,  I 

define  the  maximum  of  the  distance  vector 


s 


SRD 


=  Max[£s,„(x,,>>,)].^, 


(16) 


which  is  readily  implemented  in  FORTRAN  90  using  an  elemental  function. 

SRD=Max(SymRelDif((x,y)) 

Where  the  SymRelDif  function  is 
Elemental  Function  SymRelDif(x,y) 

Real::Intent(in)::  x,y 
ReahrSymRelDif 
If  (x=y)  then 
SymRelDif=0 
Else 

SymRelDif  =Abs(x-y)/((abs(x)+abs(y))/2) 

End  if 

End  Function 

Implementing  the  SymRelDif  Fnnction  in  Fortran 
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